
library(tidyverse)
library(haven)
library(fixest)
library(ggsci)
library(ggthemes)
library(cowplot)
library(MetBrewer)
library(PNWColors)


# Data -----

ecuador_df <- read_rds("~/Downloads/ecuador_inpraise_df.rds")
mortality_long <- read_rds("~/Downloads/ecuador_subsidy_boots_mortality_long.rds")



india_kg_pm_scenario_1_yr <- read_rds("~/Downloads/india_kg_pm_scenarios_long.rds")
india_density_2019_2023 <- read_rds("~/Downloads/india_subsidies_kg_density_2019_2023.rds")

kenya_kg_lpg <- read_rds("~/Downloads/kenya_kg_pm_scenario_1_yr_long.rds")
kenya_density_2019_2023 <- read_rds("~/Downloads/kenya_kg_vat_density_2019_2023.rds")


# Main figure -------

# ECUADOR -------

# a. Clean fuel scale up ---------

ecuador_trajectories_fig <- 
  ggplot() + 
  geom_line(data=ecuador_df, aes(x=Year, y=cf), color="black", size=0.9) +
  # geom_line(data=cf_scenario, aes(x=Year, y=cf_ctr10yr), color="dodgerblue4") + 
  geom_line(data=ecuador_df, aes(x=Year, y=cf_ctr20yr), color="black", linetype="dotted", size=0.9) +
  annotate("text", x=1990, y=0.75, label="Observed") + 
  # annotate("text", x=1991, y=.45, label="10 y behind", color="dodgerblue4") + 
  annotate("text", x=2006, y=0.3, label="Absent LPG subsidy", color="black") + 
  scale_y_continuous(labels=scales::percent_format(),
                     breaks=c(0, .25, .5, .75, 1)) + 
  xlab("") + ylab("Primary clean cooking fuel use") + 
  coord_cartesian(ylim=c(0,1), clip="off") + 
  ggtitle("Counterfactual clean cooking fuel use in Ecuador") + 
  theme_linedraw() +
  theme(
    axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    axis.title.x = element_blank(),
  )


# b. Air pollution exposure contrasts --------

ecuador_pm_contrasts <- 
  mortality_long %>% 
  dplyr::select(pm_polluting, pm_clean) %>% 
  pivot_longer(everything()) %>% 
  filter(!is.na(value))


label_ug_m3 <- function(x) {
  paste0(x, " µg/m³")
}

ecaudor_pm_fig <- 
  ggplot(ecuador_pm_contrasts, 
       aes(x=name, y=value, color=name, fill=name)) + 
  # geom_jitter(alpha=0.1, width=0.25, size=0.25) +
  geom_boxplot(width=0.15, color="black",outlier.shape = NA) + 
  stat_summary(fun = "mean", colour = "white", size = 2, geom = "point") + 
  scale_fill_manual(values=c(met.brewer("VanGogh1", 7)[6], met.brewer("VanGogh1")[2])) +
  scale_color_manual(values=c(met.brewer("VanGogh1", 7)[6], met.brewer("VanGogh1")[2])) +
  # scale_color_aaas() +
  coord_cartesian(clip="off") + 
  scale_y_continuous(breaks=c(seq(0, 120, 25)), labels = label_ug_m3) + 
  scale_x_discrete(labels=c("Primary clean", "Primary polluting")) + 
  ggtitle("PM2.5 exposures by cooking fuel type") + 
  ylab("") + 
  theme_linedraw() +
  theme(
    axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    axis.title.x = element_blank(),
  )

  
  
# c. Changes in monetized mortality, CO2e, and costs --------

ecuador_summary_mortality <- 
  mortality_long %>% 
  filter(!is.na(mortality_change_20yr)) %>% 
  group_by(model) %>% 
  mutate(mortality_cumulative = cumsum(mortality_change_20yr),
         mortality_change_VSL_20yr_cumulative = cumsum(mortality_change_VSL_20yr)) %>% 
  group_by(Year) %>%
  summarize(VSL_mean = mean(mortality_change_VSL_20yr, na.rm=T),
            deaths_mean = mean(mortality_change_20yr, na.rm=T),
            deaths_cumulative_mean = mean(mortality_cumulative, na.rm=T),
            VSL_cumulative_mean = mean(mortality_change_VSL_20yr_cumulative, na.rm=T)
  ) 

ecuador_summary_co2e <- 
  ecuador_co2e %>% 
  group_by(Year) %>% 
  summarize(
    co2e_mean = mean(co2e_difference, na.rm=T),
    co2e_usd_mean  = mean(co2e_usd, na.rm=T)
  ) %>% 
  mutate(
    co2e_usd_cumulative = cumsum(co2e_usd_mean)
  )

ecuador_summary_costs <- 
  read_xlsx("~/Downloads/LPG_subsidy_analysis.xlsx", skip=4) %>% 
  rename(Year = YEAR,
         cost_eia = `diff, EIA, 2023 dollars`,
         cost_fred = `diff, FRED, 2023 dollars`
  ) %>% 
  filter(!is.na(cost_eia)) %>% 
  mutate(
    cost_eia_cumsum = cumsum(as.numeric(cost_eia)*1e6),
    cost_fred_cumsum = cumsum(as.numeric(cost_fred)*1e6),
  )


ecuador_summary_df <- 
  ecuador_summary_mortality %>% 
  left_join(
    ecuador_summary_co2e
  ) %>% 
  left_join(
    ecuador_summary_costs
  ) %>% 
  mutate(net_benefits = (-VSL_cumulative_mean*1000) + 
           (-co2e_usd_cumulative*1000) + 
           -cost_fred_cumsum
           )

ecaudor_costs_benefits_fig <- 
  ggplot() + 
  geom_hline(yintercept = 0, color="black", linetype="dotted", size=0.5) + 
  # geom_line(
  #   data=ecuador_summary_df,
  #   aes(x=Year, y=net_benefits), linewidth=0.9
  # ) + 
  geom_line(
    data=ecuador_summary_df,
    aes(x=Year, y=-VSL_cumulative_mean/1e6), color=met.brewer("Ingres", 3)[1],
  ) + 
  geom_line(
    data=ecuador_summary_df,
    aes(x=Year, y=-co2e_usd_cumulative/1e6), color=met.brewer("Ingres", 3)[2], linetype="dotted"
  ) + 
  geom_line(
    data=ecuador_summary_df,
    aes(x=Year, y=-cost_fred_cumsum/1e9), color=met.brewer("Ingres", 3)[3], linetype="dashed"
  ) + 
  # annotate("text", x=2015, y=5.5e10, label="Net benefit") + 
  annotate("text", x=2015, y=40, label="Mortality benefit", color=met.brewer("Ingres", 3)[1]) + 
  annotate("text", x=2015, y=19, label="CO2e benefit", color=met.brewer("Ingres", 3)[2]) + 
  annotate("text", x=2015, y=-7, label="Subsidy cost", color=met.brewer("Ingres", 3)[3]) + 
  scale_y_continuous(labels=scales::dollar_format(suffix="B")) + 
  coord_cartesian(clip="off") + 
  ggtitle("Cumulative costs and benefits relative to no LPG subsidy") + 
  ylab("") + 
  theme_linedraw() +
  theme(
    axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    axis.title.x = element_blank(),
  )

ecaudor_costs_benefits_fig


# combine ------

combined_ecuador_figs <-
  cowplot::plot_grid(
    ecaudor_pm_fig,
    # ecuador_trajectories_fig,
    ecaudor_costs_benefits_fig,
    align="hv",
    labels=c("a", "b"),
    nrow=1
    )

combined_ecuador_figs


# INDIA -------

# a. Price elasticities and air pollution exposure contrasts --------

# KG LPG <> PM2.5 ------

india_kg_lpg_pm25 <- read_rds("~/Downloads/india_kg_lpg_pm25_long.rds") %>% 
  group_by(kg_lpg) %>% 
  summarize(
    PM25_1.mean=mean(PM25_1, na.rm=T),
    PM25_1.q025 = quantile(PM25_1, probs=.025, na.rm=T),
    PM25_1.q975 = quantile(PM25_1, probs=.975, na.rm=T),
  )


india_kg_lpg_pm_fig <- ggplot() + 
  geom_line(
    data=india_kg_lpg_pm25, 
    aes(x=kg_lpg, y=PM25_1.mean)
  ) +
  geom_ribbon(
    data=india_kg_lpg_pm25, 
    aes(x=kg_lpg, ymin=PM25_1.q025, ymax=PM25_1.q975), alpha=.25
  ) + 
  scale_y_continuous(breaks=c(seq(0, 250, 25)), labels = label_ug_m3) + 
  ggtitle("PM2.5 exposures and LPG use") + 
  annotate("text", x=115, y=84, label="Mean exposure per scenario",hjust=0.5) + 
  annotate("text", x=115, y=63.21639, label="63 ug/m3", color=met.brewer("VanGogh1",7)[7],hjust=0.5) + 
  annotate("text", x=115, y=66.68576, label="67 ug/m3", color=met.brewer("VanGogh1",7)[5],hjust=0.5) + 
  annotate("text", x=115, y=71.49596, label="72 ug/m3", color=met.brewer("VanGogh1",7)[4],hjust=0.5) + 
  annotate("text", x=115, y=76.15775, label="76 ug/m3", color=met.brewer("VanGogh1",7)[2],hjust=0.5) + 
  geom_segment(
    aes(
      y= mean(india_kg_pm_scenario_1_yr$PM25_1_1100, na.rm=T),
      yend=mean(india_kg_pm_scenario_1_yr$PM25_1_1100, na.rm=T),
      x=1,
      xend=58),
    color=met.brewer("VanGogh1",7)[2]) + 
  geom_segment(
    aes(
      y= 0,
      yend=mean(india_kg_pm_scenario_1_yr$PM25_1_1100, na.rm=T),
      x=58,
      xend=58),
    color=met.brewer("VanGogh1",7)[2])  + 
geom_segment(
    aes(
      y= mean(india_kg_pm_scenario_1_yr$PM25_1_900, na.rm=T),
      yend=mean(india_kg_pm_scenario_1_yr$PM25_1_900, na.rm=T),
      x=1,
      xend=61.5),
    color=met.brewer("VanGogh1",7)[4]) + 
  geom_segment(
    aes(
      y= 0,
      yend=mean(india_kg_pm_scenario_1_yr$PM25_1_900, na.rm=T),
      x=61.5,
      xend=61.5),
    color=met.brewer("VanGogh1",7)[4]) +
  geom_segment(
    aes(
      y= mean(india_kg_pm_scenario_1_yr$PM25_1_700, na.rm=T),
      yend=mean(india_kg_pm_scenario_1_yr$PM25_1_700, na.rm=T),
      x=1,
      xend=66.1),
    color=met.brewer("VanGogh1",7)[5]) + 
  geom_segment(
    aes(
      y= 0,
      yend=mean(india_kg_pm_scenario_1_yr$PM25_1_700, na.rm=T),
      x=66.1,
      xend=66.1),
    color=met.brewer("VanGogh1",7)[5]) +
  
  geom_segment(
    aes(
      y= mean(india_kg_pm_scenario_1_yr$PM25_1_550, na.rm=T),
      yend=mean(india_kg_pm_scenario_1_yr$PM25_1_550, na.rm=T),
      x=1,
      xend=69.5),
    color=met.brewer("VanGogh1",7)[7]) + 
  geom_segment(
    aes(
      y= 0,
      yend=mean(india_kg_pm_scenario_1_yr$PM25_1_550, na.rm=T),
      x=69.5,
      xend=69.5),
    color=met.brewer("VanGogh1",7)[7]) +
  ylab("") + 
  scale_x_continuous(
    limits=c(0, 200),
    breaks=c(0, 14*2, 14*4, 14*6, 14*8, 14*10, 14*12, 14*14),
    labels=c("0\n0", "28\n2","56\n4","84\n6","112\n8","140\n10","168\n12", "196\n14")
  )  +
  theme_linedraw() +
  theme(
    axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
  )

india_kg_lpg_pm_fig
# india_density_2019_2023

india_kg_lpg_fig <- ggplot() + 
  geom_line(data=india_density_2019_2023,
               aes(x=kg_lpg, y=pred_2023_1100_density_values), color=met.brewer("VanGogh1",7)[2]) + 
  geom_line(data=india_density_2019_2023,
            aes(x=kg_lpg, y=pred_2023_900_density_values), color=met.brewer("VanGogh1",7)[4]) + 
  geom_line(data=india_density_2019_2023,
            aes(x=kg_lpg, y=pred_2023_700_density_values), color=met.brewer("VanGogh1",7)[5]) + 
  geom_line(data=india_density_2019_2023,
            aes(x=kg_lpg, y=pred_2023_550_density_values), color=met.brewer("VanGogh1",7)[7]) + 
  annotate("text", x=150, y=.0035, label="550 INR", color=met.brewer("VanGogh1",7)[7]) + 
  annotate("text", x=150, y=.007, label="700 INR", color=met.brewer("VanGogh1",7)[5]) + 
  annotate("text", x=150, y=.0105, label="900 INR", color=met.brewer("VanGogh1",7)[4]) + 
  annotate("text", x=150, y=.014, label="1100 INR", color=met.brewer("VanGogh1",7)[2]) + 
  scale_x_continuous(
    limits=c(0, 200),
    breaks=c(0, 14*2, 14*4, 14*6, 14*8, 14*10, 14*12, 14*14),
    labels=c("0 kg\n0 cylinders", "28\n2","56\n4","84\n6","112\n8","140\n10","168\n12", "196 kg\n14 cylinders")
  ) + 
  xlab("Yearly LPG consumption") + 
  ylab("density") + 
  # xlim(0, 200) + 
  # coord_cartesian(xlim=c(0, 200)) +
  theme_linedraw() +
  theme(
    # axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black")
  )

india_kg_lpg_fig

india_kg_lpg_combined_fig <- cowplot::plot_grid(
  india_kg_lpg_pm_fig,
  india_kg_lpg_fig,
  nrow=2,
  align="v",
  rel_heights=c(1, 0.4)
)


# c. Changes in monetized mortality, CO2e, and costs --------

# we include 520 million people
# given 5-5.5 people per household
# indicates 96 million households - 104 million households
# mean refills per PMUY household somewhere between 3-4.4 (links below)
# https://thewire.in/government/in-fy23-lpg-refills-under-ujjwala-fell-12-of-beneficiaries-didnt-get-a-single-refill-govt
# and Report of the - Comptroller and Auditor General of India https://cag.gov.in/webroot/uploads/download_audit_report/2019/Report_No_14_of_2019_Performance_Audit_of_Pradhan_Mantri_Ujjwala_Yojana_Ministry_of_Petroleum_and_Natural_Gas_0.pdf
# intuitively, would result in 100 m * (3-4.4) *14.2 kg LPG per year
# 100000000*3.7*14.2 = 5.254e+09 kg per year # chosen as central estimate
# 100000000*4.4*14.2 = 6.248e+09 kg per year


india_kg_pm_scenario_1_yr_long %>% 
  group_by(Year) %>% 
  summarize(
    pred_550_kg = mean(pred_550_kg, na.rm=T),
    pred_700_kg = mean(pred_700_kg, na.rm=T),
    pred_900_kg = mean(pred_900_kg, na.rm=T),
    pred_1100_kg = mean(pred_1100_kg, na.rm=T)
  )

# our estimates are pretty close. no correction needed. 


# subsidized costs -----

india_subsidy_summary <- 
  india_kg_pm_scenario_1_yr_long %>%
  mutate(
  ) %>% 
  mutate(
    diff_900 = pred_1100_kg - pred_900_kg,
    diff_700 = pred_1100_kg - pred_700_kg,
    diff_550 = pred_1100_kg - pred_550_kg
  ) %>% 
  mutate(
    subidy_900_usd_total = (.93-.76) * diff_900,
    subidy_700_usd_total = (.93-.59) * diff_900,
    subidy_550_usd_total = (.93-.47) * diff_550
  ) %>% 
  mutate(
    subidy_900_usd_total_disc = subidy_900_usd_total / (1 + 0.09)^(Year-2023),
    subidy_700_usd_total_disc = subidy_700_usd_total / (1 + 0.09)^(Year-2023),
    subidy_550_usd_total_disc = subidy_550_usd_total / (1 + 0.09)^(Year-2023)
  ) %>% 
  group_by(model) %>% 
  mutate(
    subidy_900_usd_total_disc_cumsum = cumsum(subidy_900_usd_total_disc),
    subidy_700_usd_total_disc_cumsum = cumsum(subidy_700_usd_total_disc),
    subidy_550_usd_total_disc_cumsum = cumsum(subidy_550_usd_total_disc)
  ) %>% 
  group_by(Year) %>% 
  summarise( 
    subidy_900_usd = mean(subidy_900_usd_total_disc_cumsum, na.rm=T),
    subidy_700_usd = mean(subidy_700_usd_total_disc_cumsum, na.rm=T),
    subidy_550_usd = mean(subidy_550_usd_total_disc_cumsum, na.rm=T)
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 8, 10)
  ) 

# mortality benefits -----

# a similar, but different correction must be applied to mortality data
# we predict between 1.1-1.2 million premature deaths due to HAP among PMUY beneficiares
# GBD estimates 600k deaths due to HAP in India overall.
# we consider that PMUY beneficiares likely comprise 2/3 - 5/6 of these deaths
# we apply a correction:

# 1107205 / (606889*(2/3)) = 2.736592
# 1107205 / (606889*(5/6)) = 2.189273

india_kg_pm_scenario_1_yr_long <- india_kg_pm_scenario_1_yr

india_kg_pm_scenario_1_yr_long %>% 
  mutate(
    pred_900_deaths_difference = pred_1100_deaths1 - pred_900_deaths1,
    pred_700_deaths_difference = pred_1100_deaths1 - pred_700_deaths1,
    pred_550_deaths_difference = pred_1100_deaths1 - pred_550_deaths1
  ) %>% 
  group_by(model) %>% 
  mutate(
    pred_900_deaths_difference_cumsum = cumsum(pred_900_deaths_difference),
    pred_700_deaths_difference_cumsum = cumsum(pred_700_deaths_difference),
    pred_550_deaths_difference_cumsum = cumsum(pred_550_deaths_difference)
  ) %>% 
  group_by(
    Year
  ) %>% 
  summarize(
    pred_1100_deaths1 = mean(pred_1100_deaths1),
    pred_550_deaths1 = mean(pred_550_deaths1),
    pred_900_deaths_difference = mean(pred_900_deaths_difference, na.rm=T),
    pred_700_deaths_difference = mean(pred_700_deaths_difference, na.rm=T),
    pred_550_deaths_difference = mean(pred_550_deaths_difference, na.rm=T),
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 6, 8)
  )


india_mortality_summary <- 
  india_kg_pm_scenario_1_yr_long %>% 
  # filter(!is.na(kg_lpg_1100)) %>% 
  group_by(model) %>% 
  mutate(
    # pred_1100_deaths_difference_VSL1 = cumsum(pred_1100_deaths_difference_VSL1),
    pred_900_deaths_difference_VSL1 = cumsum(pred_900_deaths_difference_VSL1_discounted/2.736592),
    pred_700_deaths_difference_VSL1 = cumsum(pred_700_deaths_difference_VSL1_discounted/2.736592),
    pred_550_deaths_difference_VSL1 = cumsum(pred_550_deaths_difference_VSL1_discounted/2.736592)
  ) %>% 
  group_by(
    Year
  ) %>% 
  summarize(
    pred_900_deaths_difference_VSL1 = mean(pred_900_deaths_difference_VSL1, na.rm=T),
    pred_700_deaths_difference_VSL1 = mean(pred_700_deaths_difference_VSL1, na.rm=T),
    pred_550_deaths_difference_VSL1 = mean(pred_550_deaths_difference_VSL1, na.rm=T),
    # kg_lpg_550_cumsum = mean(kg_lpg_550_cumsum, na.rm=T),
    # 
    # kg_lpg_900_diff = mean(kg_lpg_1100_cumsum-kg_lpg_900_cumsum, na.rm=T),
    # kg_lpg_700_diff = mean(kg_lpg_1100_cumsum-kg_lpg_700_cumsum, na.rm=T),
    # kg_lpg_550_diff = mean(kg_lpg_1100_cumsum-kg_lpg_550_cumsum, na.rm=T)
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 6, 8)
  )


# co2e benefits -------


india_co2e <- 
  india_kg_pm_scenario_1_yr_long %>% 
  filter(!is.na(pred_1100_kg)) %>% 
  mutate(
    lpg_1100_mj = pred_1100_kg *  45, 
    lpg_900_mj = pred_900_kg * 45, 
    lpg_700_mj = pred_700_kg * 45, 
    lpg_550_mj = pred_550_kg * 45, 
  ) %>% 
  mutate(
    lpg_1100_mj_deliv = lpg_1100_mj * .5,
    lpg_900_mj_deliv = lpg_900_mj * .5,
    lpg_700_mj_deliv = lpg_700_mj * .5,
    lpg_550_mj_deliv = lpg_550_mj * .5
  ) %>% 
  mutate(
    delta_lpg_900_mj_deliv = lpg_1100_mj_deliv - lpg_900_mj_deliv,
    delta_lpg_700_mj_deliv = lpg_1100_mj_deliv - lpg_700_mj_deliv,
    delta_lpg_550_mj_deliv = lpg_1100_mj_deliv - lpg_550_mj_deliv
  ) %>% 
  mutate(
    delta_biomass_mj_900_deliv = delta_lpg_900_mj_deliv * (.5/.15),
    delta_biomass_mj_700_deliv = delta_lpg_700_mj_deliv * (.5/.15),
    delta_biomass_mj_550_deliv = delta_lpg_550_mj_deliv * (.5/.15)
  ) %>% 
  mutate(
    delta_lpg_co2e_900 = delta_lpg_900_mj_deliv * 173,
    delta_biomass_co2e_900 = delta_biomass_mj_900_deliv * 275,
    
    delta_lpg_co2e_700 = delta_lpg_700_mj_deliv * 173,
    delta_biomass_co2e_700 = delta_biomass_mj_700_deliv * 275,
    
    delta_lpg_co2e_550 = delta_lpg_550_mj_deliv * 173,
    delta_biomass_co2e_550 = delta_biomass_mj_550_deliv * 275
  )  %>% 
  mutate(
    delta_co2e_900 = (delta_lpg_co2e_900 + delta_biomass_co2e_900)/ 1000000000000,# grams to million kg (megatonne)
    delta_co2e_700 = (delta_lpg_co2e_700 + delta_biomass_co2e_700)/ 1000000000000,
    delta_co2e_550 = (delta_lpg_co2e_550 + delta_biomass_co2e_550)/ 1000000000000
  )%>% 
  mutate(
    delta_co2e_900_usd = delta_co2e_900 * 203, 
    delta_co2e_700_usd = delta_co2e_700 * 203,
    delta_co2e_550_usd = delta_co2e_550 * 203
  ) %>% 
  mutate(
    delta_co2e_900_usd_discounted = delta_co2e_900_usd / (1 + 0.09)^(Year-2020), # 203 usd per tone in 2020, discount to 2020
    delta_co2e_700_usd_discounted = delta_co2e_700_usd / (1 + 0.09)^(Year-2020),
    delta_co2e_550_usd_discounted = delta_co2e_550_usd / (1 + 0.05)^(Year-2020)
  ) 

india_co2e_summary <- 
  india_co2e %>% 
  group_by(model) %>% 
  mutate(
    delta_co2e_900_usd_discounted_cumsum = cumsum(delta_co2e_900_usd_discounted),
    delta_co2e_700_usd_discounted_cumsum = cumsum(delta_co2e_700_usd_discounted),
    delta_co2e_550_usd_discounted_cumsum = cumsum(delta_co2e_550_usd_discounted)
  ) %>% 
  group_by(
    Year
  ) %>% 
  summarize(
    delta_co2e_900_usd_discounted_cumsum = mean(delta_co2e_900_usd_discounted_cumsum),
    delta_co2e_700_usd_discounted_cumsum = mean(delta_co2e_700_usd_discounted_cumsum),
    delta_co2e_550_usd_discounted_cumsum = mean(delta_co2e_550_usd_discounted_cumsum)
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 12, 14)
  )


# combined costs and benefits ------
india_costs_benefits_df <-
  india_co2e_summary %>% 
  rename(co2e=value) %>% 
  left_join(
    india_mortality_summary %>% 
      rename(deaths=value), by=c("Year", "price")
  ) %>% 
  left_join(
    india_subsidy_summary %>% 
      rename(subsidy_cost=value), by=c("Year", "price")
  )


india_cost_benefit_summary_fig <- 
  ggplot() + 
  geom_line(data=india_costs_benefits_df, aes(x=Year, y=-co2e/1000, color=price, group=price), linetype="dotted") + # scale to billions
  geom_line(data=india_costs_benefits_df, aes(x=Year, y=subsidy_cost/1000000000, color=price, group=price), linetype="dashed") +# scale to billions
  geom_line(data=india_costs_benefits_df, aes(x=Year, y=deaths/1000000000, color=price, group=price)) +# scale to billions
  ggtitle("Cumulative costs and benefits relative to market cost (1100 INR)") + 
  geom_hline(yintercept=0, linesize=0.25, color="black") + 
  scale_color_manual(values = c(met.brewer("VanGogh1",7)[7],
                                met.brewer("VanGogh1",7)[5],
                                met.brewer("VanGogh1",7)[4])) + 
  annotate("text", x=2026, y=150, label="550 INR", color=met.brewer("VanGogh1",7)[7]) + 
  annotate("text", x=2026, y=165, label="700 INR", color=met.brewer("VanGogh1",7)[5]) + 
  annotate("text", x=2026, y=180, label="900 INR", color=met.brewer("VanGogh1",7)[4]) + 
  theme_linedraw() +
  geom_segment(aes(x=2023,xend=2023.5,y=180,yend=180)) + 
  annotate("text", x=2024, y=180, label="Mortality", color="black",hjust=0.2) +
  geom_segment(aes(x=2023,xend=2023.5,y=165,yend=165), linetype="dotted") + 
  annotate("text", x=2024, y=165, label="CO2e", color="black",hjust=0.2) +
  geom_segment(aes(x=2023,xend=2023.5,y=150,yend=150), linetype="dashed") + 
  annotate("text", x=2024, y=150, label="Subsidy", color="black",hjust=0.2) +
  # coord_cartesian(xlim=c(2023, 2031), clip = "off") + 
  scale_y_continuous(labels = scales::dollar_format(suffix = "B")) + 
  theme_linedraw() +
  theme(
    axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    legend.title = element_blank(),
    axis.title.x = element_blank(),
  )

india_cost_benefit_summary_fig

# d. Summarized net benefits ----------

combined_india_figs <- 
  plot_grid(
  india_kg_lpg_combined_fig,
    # india_kg_summary_fig,
    india_cost_benefit_summary_fig,
    nrow=1,
  labels=c("c", "d")
  )

combined_india_figs

# KENYA -------

kenya_kg_pm_scenario_1_yr_long <- read_rds("~/Downloads/kenya_kg_pm_scenario_1_yr_long.rds")


# a. Price elasticities and air pollution exposure contrasts --------

kenya_kg_lpg_pm25 <- read_rds("~/Downloads/kenya_kg_lpg_pm25_long.rds") %>% 
  group_by(kg_lpg) %>% 
  summarize(
    PM25_1.mean=mean(PM25_1, na.rm=T),
    PM25_1.q025 = quantile(PM25_1, probs=.025, na.rm=T),
    PM25_1.q975 = quantile(PM25_1, probs=.975, na.rm=T),
  )


kenya_lpg_pm25 <- 
  ggplot() +
  geom_line(data=kenya_kg_lpg_pm25, aes(x=kg_lpg, y=PM25_1.mean)) +
  geom_ribbon(data=kenya_kg_lpg_pm25, aes(x=kg_lpg, ymin=PM25_1.q025, ymax=PM25_1.q975), alpha=0.25) +
  scale_y_continuous(breaks=c(seq(0, 300, 25)), labels = label_ug_m3) + 
  coord_cartesian(ylim=c(0, 300), xlim=c(0, 55)) +
  
  annotate("text", x=30, y=200, label="Mean exposure per scenario", hjust=0.5) + 
  annotate("text", x=30, y=187, label="186 ug/m3", color=met.brewer("VanGogh1",4)[1], hjust=0.5) +
  annotate("text", x=30, y=172, label="177 ug/m3", color=met.brewer("VanGogh1",7)[5], hjust=0.5) + 
  annotate("text", x=30, y=160, label="160 ug/m3", color=met.brewer("VanGogh1",7)[7], hjust=0.5) +  
  geom_segment(
    aes(
      y= 0,
      yend=mean(kenya_kg_pm_scenario_1_yr_long$PM25_1_16, na.rm=T),
      x=13.1,
      xend=13.1),
    color=met.brewer("VanGogh1",4)[1]) +
  geom_segment(
    aes(
      y= mean(kenya_kg_pm_scenario_1_yr_long$PM25_1_16, na.rm=T),
      yend=mean(kenya_kg_pm_scenario_1_yr_long$PM25_1_16, na.rm=T),
      x=1,
      xend=13.1),
    color=met.brewer("VanGogh1",4)[1]) + 
  
  geom_segment(
    aes(
      y= mean(kenya_kg_pm_scenario_1_yr_long$PM25_1_8, na.rm=T),
      yend=mean(kenya_kg_pm_scenario_1_yr_long$PM25_1_8, na.rm=T),
      x=1,
      xend=14.1),
    color=met.brewer("VanGogh1",7)[5]) + 
  geom_segment(
    aes(
      y= 0,
      yend=mean(kenya_kg_pm_scenario_1_yr_long$PM25_1_8, na.rm=T),
      x=14.1,
      xend=14.1),
    color=met.brewer("VanGogh1",7)[5]) +
  
  geom_segment(
    aes(
      y= 0,
      yend=mean(kenya_kg_pm_scenario_1_yr_long$PM25_1_observed, na.rm=T),
      x=15.99,
      xend=15.99),
    color=met.brewer("VanGogh1",7)[7]) +
  
  geom_segment(
    aes(
      y= mean(kenya_kg_pm_scenario_1_yr_long$PM25_1_observed, na.rm=T),
      yend=mean(kenya_kg_pm_scenario_1_yr_long$PM25_1_observed, na.rm=T),
      x=1,
      xend=15.99),
    color=met.brewer("VanGogh1",7)[7]) + 
  
  
  ggtitle("PM2.5 exposure and LPG use") + 
  ylab("Mean PM2.5 exposure ug/m3") + 
  xlab("LPG consumption kg/capita/year") + 
  theme_linedraw() +
  theme(
    axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
  )

kenya_lpg_pm25

kenya_density_2019_2023


kenya_lpg_densities_fig <- 
  ggplot() + 
  geom_line(data=kenya_density_2019_2023,
               aes(x=kg_lpg, y=observed_density_values_1), color=met.brewer("VanGogh1",7)[7]) + 
  geom_line(data=kenya_density_2019_2023,
            aes(x=kg_lpg, y=pred_2023_1_16_density_values_1), color=met.brewer("VanGogh1",7)[5]) + 
  geom_line(data=kenya_density_2019_2023,
            aes(x=kg_lpg, y=pred_2023_1_8_density_values_1), color=met.brewer("VanGogh1",4)[1]) + 
            
  annotate("text", x=40, y=.015, label="No VAT", color=met.brewer("VanGogh1",7)[7]) + 
  annotate("text", x=40, y=.03, label="8% VAT", color=met.brewer("VanGogh1",7)[5]) + 
  annotate("text", x=40, y=.045, label="16% VAT", color=met.brewer("VanGogh1",4)[1]) + 
  scale_x_continuous(
    limits=c(0, 55)
  ) +
  xlab("LPG consumption kg/capita/year") + 
  ylab("density") + 
  coord_cartesian(xlim=c(0, 55)) + 
  theme_linedraw() +
  theme(
    # axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    # axis.text.x = element_blank(),
    # axis.title.x = element_blank(),
  )


kenya_kg_lpg_combined_fig <- cowplot::plot_grid(
  kenya_lpg_pm25,
  kenya_lpg_densities_fig,
  nrow=2,
  align="v",
  rel_heights=c(1, 0.4)
)

kenya_kg_lpg_combined_fig

kenya_pm25_distribution_fig <- 
  ggplot() + 
  geom_boxplot(data=kenya_kg_pm_scenario_1_yr_long,
            aes(x="No VAT", y=PM25_1_observed), fill=met.brewer("VanGogh1",7)[7], width=0.15, color="black",outlier.shape = NA) + 
  geom_boxplot(data=kenya_kg_pm_scenario_1_yr_long,
            aes(x="8% VAT", y=PM25_1_8), fill=met.brewer("VanGogh1",7)[5], width=0.15, color="black",outlier.shape = NA) + 
  geom_boxplot(data=kenya_kg_pm_scenario_1_yr_long,
            aes(x="16% VAT", y=PM25_1_16), fill=met.brewer("VanGogh1",4)[1], width=0.15, color="black",outlier.shape = NA) + 
  
  stat_summary(data=kenya_kg_pm_scenario_1_yr_long,
               aes(x="No VAT", y=PM25_1_observed), fun = "mean", colour = "white", size = 2, geom = "point") + 
  stat_summary(data=kenya_kg_pm_scenario_1_yr_long,
               aes(x="8% VAT", y=PM25_1_8), fun = "mean", colour = "white", size = 2, geom = "point") + 
  stat_summary(data=kenya_kg_pm_scenario_1_yr_long,
               aes(x="16% VAT", y=PM25_1_16), fun = "mean", colour = "white", size = 2, geom = "point") + 
  
  theme_linedraw() +
  theme(
    axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    # axis.text.x = element_blank(),
    axis.title.x = element_blank(),
  )


# kg use scenarios -------

kenya_kg_summary <- 
  kenya_kg_pm_scenario_1_yr_long %>% 
  group_by(model) %>% 
  mutate(
    kg_lpg_consumed_observed_cumsum = cumsum(kg_lpg_consumed_observed),
    kg_lpg_consumed_1_16_cumsum = cumsum(kg_lpg_consumed_1_16),
    kg_lpg_consumed_1_8_cumsum = cumsum(kg_lpg_consumed_1_8)
  ) %>% 
  ungroup() %>% 
  group_by(
    Year
  ) %>% 
  summarize(
    kg_lpg_consumed_observed_cumsum = mean(kg_lpg_consumed_observed_cumsum, na.rm=T),
    kg_lpg_consumed_1_16_cumsum = mean(kg_lpg_consumed_1_16_cumsum, na.rm=T),
    kg_lpg_consumed_1_8_cumsum = mean(kg_lpg_consumed_1_8_cumsum, na.rm=T)
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = factor(name, levels=c("kg_lpg_consumed_observed_cumsum", 
                                  "kg_lpg_consumed_1_8_cumsum",
                                  "kg_lpg_consumed_1_16_cumsum"))
  )

kenya_kg_summary_fig <- 
  ggplot(
    kenya_kg_summary, aes(x=Year, y=value, color=price, group=price)
  ) + 
  geom_line() +
  ggtitle("Cumulative LPG consumption") + 
  scale_color_manual(values = c(met.brewer("VanGogh1",7)[7],
                                met.brewer("VanGogh1",7)[5],
                                met.brewer("VanGogh1",4)[1])) + 
  annotate("text", x=2026, y=5e9, label="No VAT", color=met.brewer("VanGogh1",7)[7]) + 
  annotate("text", x=2026, y=4.5e9, label="8% VAT", color=met.brewer("VanGogh1",7)[5]) + 
  annotate("text", x=2026, y=4e9, label="16% VAT", color=met.brewer("VanGogh1",4)[1]) + 
  theme_linedraw() +
  theme(
    axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    axis.title.x = element_blank(),
  )


kenya_kg_summary_fig

# c. Changes in monetized mortality, CO2e, and costs --------

# subsidized costs -----


base_price = 1.49138 # current price in kenya

kenya_subsidy_summary <- 
  kenya_kg_pm_scenario_1_yr_long %>%
  mutate(
    kg_lpg_consumed_observed =  kg_lpg_consumed_observed / 1.6,
    kg_lpg_consumed_1_16 =  kg_lpg_consumed_1_16 / 1.6,
    kg_lpg_consumed_1_8 =  kg_lpg_consumed_1_8 / 1.6
  ) %>% 
  mutate(
    diff_1_16 = kg_lpg_consumed_observed - kg_lpg_consumed_1_16,
    diff_1_8 = kg_lpg_consumed_observed - kg_lpg_consumed_1_8
  ) %>% 
  mutate(
    subidy_1_16_usd_total = (base_price*.16 -base_price) * diff_1_16,
    subidy_1_8_usd_total = (base_price*.08 -base_price) * diff_1_8
  ) %>% 
  mutate(
    subidy_1_16_usd_total_disc = subidy_1_16_usd_total / (1 + 0.09)^(Year-2023),
    subidy_1_8_usd_total_disc = subidy_1_8_usd_total / (1 + 0.09)^(Year-2023)
  ) %>% 
  group_by(model) %>% 
  mutate(
    subidy_1_16_usd_total_disc_cumsum = cumsum(subidy_1_16_usd_total_disc),
    subidy_1_8_usd_total_disc_cumsum = cumsum(subidy_1_8_usd_total_disc)
  ) %>% 
  group_by(Year) %>% 
  summarise( 
    subidy_1_16_usd = mean(subidy_1_16_usd_total_disc_cumsum, na.rm=T),
    subidy_1_8_usd = mean(subidy_1_8_usd_total_disc_cumsum, na.rm=T)
  ) 

# mortality benefits -----
colnames(kenya_kg_pm_scenario_1_yr_long)

kenya_mortality_summary <- 
  kenya_kg_pm_scenario_1_yr_long %>% 
  # filter(!is.na(kg_lpg_1100)) %>% 
  group_by(model) %>% 
  mutate(
    pred_1_16_deaths_difference1_VSL_discounted_cumsum = cumsum(pred_1_16_deaths_difference1_VSL_discounted),
    pred_1_8_deaths_difference1_VSL_discounted_cumsum = cumsum(pred_1_8_deaths_difference1_VSL_discounted),
  ) %>% 
  group_by(
    Year
  ) %>% 
  summarize(
    pred_1_16_deaths_difference1_VSL_discounted_cumsum = mean(pred_1_16_deaths_difference1_VSL_discounted_cumsum, na.rm=T),
    pred_1_8_deaths_difference1_VSL_discounted_cumsum = mean(pred_1_8_deaths_difference1_VSL_discounted_cumsum, na.rm=T)
  ) 

# co2e benefits -------

kenya_co2e <- 
  kenya_kg_pm_scenario_1_yr_long %>%  
  mutate(
    kg_lpg_consumed_observed_mj = kg_lpg_consumed_observed / 1.6 * 45,
    kg_lpg_consumed_1_16_mj = kg_lpg_consumed_1_16  / 1.6* 45,
    kg_lpg_consumed_1_8_mj = kg_lpg_consumed_1_8 / 1.6 * 45
  ) %>% 
  mutate(
    delta_lpg_mj_1_16 = kg_lpg_consumed_1_16_mj - kg_lpg_consumed_observed_mj,
    delta_lpg_mj_1_8 = kg_lpg_consumed_1_8_mj - kg_lpg_consumed_observed_mj
  ) %>% 
  mutate(
    biomass_mj_1_16 = delta_lpg_mj_1_16 * (.5/.15),
    biomass_mj_1_8 = delta_lpg_mj_1_8 * (.5/.15)
  ) %>% 
  mutate(
    delta_lpg_co2e_1_16 = delta_lpg_mj_1_16*173,
    biomass_co2e_1_16 = biomass_mj_1_16*395,
    
    delta_lpg_co2e_1_8 = delta_lpg_mj_1_8*173,
    biomass_co2e_1_8 = biomass_mj_1_8*395
  )  %>% 
  mutate(
    delta_co2e_1_16 = (delta_lpg_co2e_1_16 + biomass_co2e_1_16)/ 1000000000000,
    delta_co2e_1_8 = (delta_lpg_co2e_1_8 + biomass_co2e_1_8)/ 1000000000000
  )%>% 
  mutate(
    delta_co2e_1_16_usd = delta_co2e_1_16 * 203 * 1e6, # from Burke et al. in 2020 203 USD / ton
    delta_co2e_1_8_usd = delta_co2e_1_8 * 203 * 1e6
  ) %>% 
  mutate(
    delta_co2e_1_16_usd_discounted = delta_co2e_1_16_usd / (1 + 0.05)^(Year-2020),
    delta_co2e_1_8_usd_discounted = delta_co2e_1_8_usd / (1 + 0.05)^(Year-2020)
  ) %>% 
  group_by(model) %>% 
  mutate(
    delta_co2e_1_16_usd_discounted_cumsum = cumsum(delta_co2e_1_16_usd_discounted),
    delta_co2e_1_8_usd_discounted_cumsum = cumsum(delta_co2e_1_8_usd_discounted)
  )

kenya_co2e_summary <- 
  kenya_co2e %>% 
  group_by(
    Year
  ) %>% 
  summarize(
    delta_co2e_1_16.mean = mean(delta_co2e_1_16, na.rm=T),
    delta_co2e_1_16.q025 = quantile(delta_co2e_1_16, probs=.025, na.rm=T),
    delta_co2e_1_16.q975 = quantile(delta_co2e_1_16, probs=.975, na.rm=T),
    delta_co2e_1_8 = mean(delta_co2e_1_8, na.rm=T),
    
    delta_co2e_1_16_usd_discounted = mean(delta_co2e_1_16_usd_discounted_cumsum, na.rm=T),
    delta_co2e_1_16_usd_discounted.q025 = quantile(delta_co2e_1_16_usd_discounted_cumsum, probs=.025, na.rm=T),
    delta_co2e_1_16_usd_discounted.q975 = quantile(delta_co2e_1_16_usd_discounted_cumsum, probs=.975, na.rm=T),
    
    delta_co2e_1_8_usd_discounted = mean(delta_co2e_1_8_usd_discounted_cumsum, na.rm=T)
  )


# combined costs and benefits ------
kenya_costs_benefits_df <-
  kenya_co2e_summary %>% 
  left_join(
    kenya_mortality_summary, by=c("Year")
  ) %>% 
  left_join(
    kenya_subsidy_summary, by=c("Year")
  ) %>% 
  mutate(
    delta_co2e_1_8_usd_discounted = delta_co2e_1_16_usd_discounted - delta_co2e_1_8_usd_discounted,
    subidy_1_8_usd = subidy_1_16_usd - subidy_1_8_usd,
    pred_1_8_deaths_difference1_VSL_discounted_cumsum = pred_1_16_deaths_difference1_VSL_discounted_cumsum - pred_1_8_deaths_difference1_VSL_discounted_cumsum
  )


kenya_cost_benefit_summary_fig <- 
  ggplot() + 
  geom_line(data=kenya_costs_benefits_df, aes(x=Year, y=-delta_co2e_1_16_usd_discounted/1e9), linetype="dotted", color=met.brewer("VanGogh1",7)[7]) +
  geom_line(data=kenya_costs_benefits_df, aes(x=Year, y=-delta_co2e_1_8_usd_discounted/1e9), linetype="dotted", color=met.brewer("VanGogh1",7)[5]) +
  
  geom_line(data=kenya_costs_benefits_df, aes(x=Year, y=subidy_1_16_usd/1e9), linetype="dashed", color=met.brewer("VanGogh1",7)[7]) +
  geom_line(data=kenya_costs_benefits_df, aes(x=Year, y=subidy_1_8_usd/1e9), linetype="dashed", color=met.brewer("VanGogh1",7)[5]) +
  
  geom_line(data=kenya_costs_benefits_df, aes(x=Year, y=-pred_1_16_deaths_difference1_VSL_discounted_cumsum/1e9), color=met.brewer("VanGogh1",7)[7]) +
  geom_line(data=kenya_costs_benefits_df, aes(x=Year, y=-pred_1_8_deaths_difference1_VSL_discounted_cumsum/1e9), color=met.brewer("VanGogh1",7)[5]) +
  
  # ylab("Yearly LPG consumption") + 
  ggtitle("Cumulative costs and benefits relative to 16% VAT") + 
  geom_hline(yintercept=0, linesize=0.25, color="black") + 
  annotate("text", x=2026, y=8.75, label="No VAT", color=met.brewer("VanGogh1",7)[7]) + 
  annotate("text", x=2026, y=9.5, label="8% VAT", color=met.brewer("VanGogh1",7)[5]) + 
  theme_linedraw() +
  geom_segment(aes(x=2023,xend=2023.5,y=9.5,yend=9.5)) + 
  annotate("text", x=2024, y=9.5, label="Mortality", color="black",hjust=0.2) +
  geom_segment(aes(x=2023,xend=2023.5,y=8.75,yend=8.75), linetype="dotted") + 
  annotate("text", x=2024, y=8.75, label="CO2e", color="black",hjust=0.2) +
  geom_segment(aes(x=2023,xend=2023.5,y=8,yend=8), linetype="dashed") + 
  annotate("text", x=2024, y=8, label="Subsidy", color="black",hjust=0.2) +
  
  # coord_cartesian(xlim=c(2023, 2031), clip = "off") + 
  scale_y_continuous(labels = scales::dollar_format(suffix = "B")) + 
  
  theme_linedraw() +
  theme(
    axis.title.y = element_blank(),
    panel.grid = element_blank(),
    # axis.text = element_text(size=10, color="black"),
    legend.position = "none",
    legend.title = element_blank(),
    axis.title.x = element_blank(),
  )

kenya_cost_benefit_summary_fig

# d. Summarized net benefits ----------

combined_kenya_figs <- 
  plot_grid(
    kenya_kg_lpg_combined_fig,
    # kenya_kg_summary_fig,
    kenya_cost_benefit_summary_fig,
    nrow=1,
    labels=c("e", "f")
  )


# COMBINE ALL FIGURES --------

combined_figs <-
  plot_grid(
    combined_ecuador_figs,
    combined_india_figs,
    combined_kenya_figs,
    nrow=3,
    align="hv"
  )

combined_figs_annotate <- 
  plot_grid(
    NULL, combined_figs,
    nrow=1, 
    rel_widths = c(0.025, 1)
  ) + 
  annotate("text", x=.015, y=0.85, label="Ecuador", angle=90, size=6,fontface="bold")  + 
  annotate("text", x=.015, y=0.525, label="India", angle=90, size=6,fontface="bold")  + 
  annotate("text", x=.015, y=0.2, label="Kenya", angle=90, size=6,fontface="bold") 


cowplot::ggsave2("~/Desktop/combined_inpraise_fig_r1.pdf",
                 combined_figs_annotate,
                 width = 300,
                 height = 350,
                 units = "mm",
                 dpi = 300)



# CALCULATE STATISTICS FOR PAPER -------


mort_df_india <- 
  india_kg_pm_scenario_1_yr_long %>% 
  mutate(
    pred_900_deaths_difference = pred_1100_deaths1 - pred_900_deaths1,
    pred_700_deaths_difference = pred_1100_deaths1 - pred_700_deaths1,
    pred_550_deaths_difference = pred_1100_deaths1 - pred_550_deaths1
  ) %>% 
  group_by(model) %>% 
  mutate(
    pred_900_deaths_difference_cumsum = cumsum(pred_900_deaths_difference),
    pred_700_deaths_difference_cumsum = cumsum(pred_700_deaths_difference),
    pred_550_deaths_difference_cumsum = cumsum(pred_550_deaths_difference)
  ) %>% 
  group_by(
    Year
  ) %>% 
  summarize(
    pred_1100_deaths1 = mean(pred_1100_deaths1),
    pred_550_deaths1 = mean(pred_550_deaths1),
    pred_900_deaths_difference_cumsum = mean(pred_900_deaths_difference_cumsum, na.rm=T),
    pred_700_deaths_difference_cumsum = mean(pred_700_deaths_difference_cumsum, na.rm=T),
    pred_550_deaths_difference_cumsum = mean(pred_550_deaths_difference_cumsum, na.rm=T),
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 6, 8)
  )

co2_df_india <- 
  india_co2e %>% 
  group_by(model) %>% 
  mutate(
    delta_co2e_900_usd_discounted_cumsum = cumsum(delta_co2e_900_usd_discounted),
    delta_co2e_700_usd_discounted_cumsum = cumsum(delta_co2e_700_usd_discounted),
    delta_co2e_550_usd_discounted_cumsum = cumsum(delta_co2e_550_usd_discounted),
    
    delta_co2e_900_cumsum = cumsum(delta_co2e_900),
    delta_co2e_700_cumsum = cumsum(delta_co2e_700),
    delta_co2e_550_cumsum = cumsum(delta_co2e_550),
  ) %>% 
  group_by(
    Year
  ) %>% 
  summarize(
    delta_co2e_900_usd_discounted_cumsum = mean(delta_co2e_900_usd_discounted_cumsum),
    delta_co2e_700_usd_discounted_cumsum = mean(delta_co2e_700_usd_discounted_cumsum),
    delta_co2e_550_usd_discounted_cumsum = mean(delta_co2e_550_usd_discounted_cumsum),
    
    delta_co2e_900_cumsum = mean(delta_co2e_900_cumsum), # millions of tons or megatonnes
    delta_co2e_700_cumsum = mean(delta_co2e_700_cumsum),
    delta_co2e_550_cumsum = mean(delta_co2e_550_cumsum)
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 12, 14)
  )

view(co2_df_india)
view(india_costs_benefits_df)
view(kenya_costs_benefits_df)
# for ref, india emits 3.9 gigatonnes CO2e each year


kenya_kg_pm_scenario_1_yr_long %>% 
  # mutate(pred_1_16_deaths_difference1)
  group_by(model) %>% 
  mutate(
    pred_1_16_deaths_difference1_cumsum = cumsum(pred_1_16_deaths_difference1),
    pred_1_8_deaths_difference1_cumsum = cumsum(pred_1_8_deaths_difference1),
  ) %>% 
  group_by(
    Year
  ) %>% 
  summarize(
    pred_1_16_deaths_difference1_cumsum = mean(pred_1_16_deaths_difference1_cumsum, na.rm=T), # relative to 0% VAT
    pred_1_8_deaths_difference1_cumsum = mean(pred_1_8_deaths_difference1_cumsum, na.rm=T) # difference 16-8 for 16 to 8% VAT
  )  # 2.6 billion for CO2

view(kenya_costs_benefits_df)







# RESULTS PLOTS SENTIVITY ---------

# Ecuador, low VSL and high VSL and low SCC and high SCC ---------

mortality_long <- read_rds("~/Downloads/ecuador_subsidy_boots_mortality_long_63.rds")

ecuador_summary_mortality <- 
  mortality_long %>% 
  filter(!is.na(mortality_change_20yr)) %>% 
  group_by(model) %>% 
  mutate(mortality_change_VSL_20yr_cumulative = cumsum(mortality_change_VSL_20yr),
         mortality_change_VSL_20yr_cumulative_low = cumsum(mortality_change_VSL_20yr_low),
         mortality_change_VSL_20yr_cumulative_high = cumsum(mortality_change_VSL_20yr_high)
  ) %>% 
  filter(Year==2019) %>% 
  ungroup() 

ecuador_summary_co2e <- 
  ecuador_co2e %>% 
  group_by(Year) %>% 
  summarize(
    co2e_mean = mean(co2e_difference, na.rm=T),
    co2e_usd_mean  = mean(co2e_usd, na.rm=T),
    co2e_usd_mean_low  = mean(co2e_usd_low, na.rm=T),
    co2e_usd_mean_high  = mean(co2e_usd_high, na.rm=T),
  ) %>% 
  mutate(
    co2e_usd_cumulative = cumsum(co2e_usd_mean),
    co2e_usd_cumulative_low = cumsum(co2e_usd_mean_low),
    co2e_usd_cumulative_high = cumsum(co2e_usd_mean_high)
  ) %>% 
  filter(Year==2019)

ecuador_summary_costs <- 
  read_xlsx("~/Downloads/LPG_subsidy_analysis.xlsx", skip=4) %>% 
  rename(Year = YEAR,
         cost_eia = `diff, EIA, 2023 dollars`,
         cost_fred = `diff, FRED, 2023 dollars`
  ) %>% 
  filter(!is.na(cost_eia)) %>% 
  mutate(
    cost_eia_cumsum = cumsum(as.numeric(cost_eia)*1e6),
    cost_fred_cumsum = cumsum(as.numeric(cost_fred)*1e6),
  ) %>% 
  filter(Year==2019)


ecuador_mortality_benefits_fig <- 
  ggplot() + 
  geom_errorbar(data=ecuador_summary_mortality %>% 
                  summarize(
                    low = quantile((-mortality_change_VSL_20yr_cumulative/1e6), probs=.025),
                    high = quantile((-mortality_change_VSL_20yr_cumulative/1e6), probs=.975)), 
                aes(x="Mortality benefits\n(preferred VSL)", ymin=low, ymax=high), linetype="dotted", color="black", width=0.15) + 
  geom_boxplot(data=ecuador_summary_mortality, 
               aes(x="Mortality benefits\n(preferred VSL)", y=(-mortality_change_VSL_20yr_cumulative/1e6)), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[1]) + 
  stat_summary(data=ecuador_summary_mortality, 
               aes(x="Mortality benefits\n(preferred VSL)", y=(-mortality_change_VSL_20yr_cumulative/1e6)), fun = "mean", colour = "white", size = 2, geom = "point") +
  
  geom_errorbar(data=ecuador_summary_mortality %>% 
                  summarize(
                    low = quantile((-mortality_change_VSL_20yr_cumulative_low/1e6), probs=.025),
                    high = quantile((-mortality_change_VSL_20yr_cumulative_low/1e6), probs=.975)), 
                aes(x="Mortality benefits\n(low VSL)", ymin=low, ymax=high), linetype="dotted", color="black", width=0.15) + 
  geom_boxplot(data=ecuador_summary_mortality, 
               aes(x="Mortality benefits\n(low VSL)", y=(-mortality_change_VSL_20yr_cumulative_low/1e6)), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[1]) +

  stat_summary(data=ecuador_summary_mortality, 
               aes(x="Mortality benefits\n(low VSL)", y=(-mortality_change_VSL_20yr_cumulative_low/1e6)), fun = "mean", colour = "white", size = 2, geom = "point") +
  
  geom_errorbar(data=ecuador_summary_mortality %>% 
                  summarize(
                    low = quantile((-mortality_change_VSL_20yr_cumulative_high/1e6), probs=.025),
                    high = quantile((-mortality_change_VSL_20yr_cumulative_high/1e6), probs=.975)), 
                aes(x="Mortality benefits\n(high VSL)", ymin=low, ymax=high), linetype="dotted", color="black", width=0.15) + 
  geom_boxplot(data=ecuador_summary_mortality, 
               aes(x="Mortality benefits\n(high VSL)", y=(-mortality_change_VSL_20yr_cumulative_high/1e6)), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[1]) + 

  stat_summary(data=ecuador_summary_mortality, 
               aes(x="Mortality benefits\n(high VSL)", y=(-mortality_change_VSL_20yr_cumulative_high/1e6)), fun = "mean", colour = "white", size = 2, geom = "point")  + 
  
  scale_x_discrete(limits=c("Mortality benefits\n(low VSL)", "Mortality benefits\n(preferred VSL)", "Mortality benefits\n(high VSL)")) +
  
  ggtitle("Mortality related benefits") + 
  xlab("") + 
  ylab("USD (billions)") + 
  
  theme_bw()



ecuador_co2_benefits_fig <- 
  ggplot() + 
  geom_col(data=ecuador_summary_co2e, 
               aes(x="CO2e benefits\n(preferred SCC)", y=(-co2e_usd_cumulative/1e6)), width=0.15, color=met.brewer("Ingres", 3)[2], fill=met.brewer("Ingres", 3)[2]) +
  geom_col(data=ecuador_summary_co2e, 
           aes(x="CO2e benefits\n(low SCC)", y=(-co2e_usd_cumulative_low/1e6)), width=0.15, color=met.brewer("Ingres", 3)[2], fill=met.brewer("Ingres", 3)[2]) +
  
  geom_col(data=ecuador_summary_co2e, 
           aes(x="CO2e benefits\n(high SCC)", y=(-co2e_usd_cumulative_high/1e6)), width=0.15, color=met.brewer("Ingres", 3)[2], fill=met.brewer("Ingres", 3)[2]) +
  scale_x_discrete(limits=c("CO2e benefits\n(low SCC)", "CO2e benefits\n(preferred SCC)", "CO2e benefits\n(high SCC)")) +
  ggtitle("CO2e related benefits") + 
  xlab("") + 
  ylab("USD (billions)") + 
  
  theme_bw()

ecuador_monetized_sensitivity_fig <- 
  cowplot::plot_grid(
  NULL, NULL,
  ecuador_mortality_benefits_fig,
  ecuador_co2_benefits_fig, nrow=2, 
  rel_heights=c(0.1, 1), 
  align="hv",
  labels="Sensitivity of estimates of monetized benefits in Ecuador to higher or lower VSLs and SCCs",
  label_x = -.5
)


# India, low VSL and high VSL and low SCC and high SCC ---------


india_mortality_summary <- 
  india_kg_pm_scenario_1_yr_long %>% 
  # filter(!is.na(kg_lpg_1100)) %>% 
  group_by(model) %>% 
  mutate(
    # pred_1100_deaths_difference_VSL1 = cumsum(pred_1100_deaths_difference_VSL1),
    pred_900_deaths_difference_VSL1 = cumsum(pred_900_deaths_difference_VSL1_discounted/2.736592),
    pred_700_deaths_difference_VSL1 = cumsum(pred_700_deaths_difference_VSL1_discounted/2.736592),
    pred_550_deaths_difference_VSL1 = cumsum(pred_550_deaths_difference_VSL1_discounted/2.736592),
    
    pred_900_deaths_difference_VSL1_low = cumsum(pred_900_deaths_difference_VSL1_low_discounted/2.736592),
    pred_700_deaths_difference_VSL1_low = cumsum(pred_700_deaths_difference_VSL1_low_discounted/2.736592),
    pred_550_deaths_difference_VSL1_low = cumsum(pred_550_deaths_difference_VSL1_low_discounted/2.736592),
    
    pred_900_deaths_difference_VSL1_high = cumsum(pred_900_deaths_difference_VSL1_high_discounted/2.736592),
    pred_700_deaths_difference_VSL1_high = cumsum(pred_700_deaths_difference_VSL1_high_discounted/2.736592),
    pred_550_deaths_difference_VSL1_high = cumsum(pred_550_deaths_difference_VSL1_high_discounted/2.736592),
    
  ) %>% 
  ungroup() %>% 
  dplyr::select(Year, 
                pred_900_deaths_difference_VSL1,
                pred_700_deaths_difference_VSL1,
                pred_550_deaths_difference_VSL1,
                
                pred_900_deaths_difference_VSL1_low,
                pred_700_deaths_difference_VSL1_low,
                pred_550_deaths_difference_VSL1_low,
                
                pred_900_deaths_difference_VSL1_high,
                pred_700_deaths_difference_VSL1_high,
                pred_550_deaths_difference_VSL1_high) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 6, 8),
    VSL = str_sub(name, nchar(name)-3,nchar(name))
  ) %>% 
  mutate(
    VSL = ifelse(VSL=="VSL1", "VSL\n(preferred)",
                 ifelse(VSL=="_low", "VSL\n(low)",
                        ifelse(VSL=="high", "VSL\n(high)",
                        NA))),
    price = paste0(price, " INR")
  ) %>% 
  filter(Year==2030)


# co2e benefits -------


india_co2e <- 
  india_kg_pm_scenario_1_yr_long %>% 
  filter(!is.na(pred_1100_kg)) %>% 
  mutate(
    lpg_1100_mj = pred_1100_kg *  45, 
    lpg_900_mj = pred_900_kg * 45, 
    lpg_700_mj = pred_700_kg * 45, 
    lpg_550_mj = pred_550_kg * 45, 
  ) %>% 
  mutate(
    lpg_1100_mj_deliv = lpg_1100_mj * .5,
    lpg_900_mj_deliv = lpg_900_mj * .5,
    lpg_700_mj_deliv = lpg_700_mj * .5,
    lpg_550_mj_deliv = lpg_550_mj * .5
  ) %>% 
  mutate(
    delta_lpg_900_mj_deliv = lpg_1100_mj_deliv - lpg_900_mj_deliv,
    delta_lpg_700_mj_deliv = lpg_1100_mj_deliv - lpg_700_mj_deliv,
    delta_lpg_550_mj_deliv = lpg_1100_mj_deliv - lpg_550_mj_deliv
  ) %>% 
  mutate(
    delta_biomass_mj_900_deliv = delta_lpg_900_mj_deliv * (.5/.15),
    delta_biomass_mj_700_deliv = delta_lpg_700_mj_deliv * (.5/.15),
    delta_biomass_mj_550_deliv = delta_lpg_550_mj_deliv * (.5/.15)
  ) %>% 
  mutate(
    delta_lpg_co2e_900 = delta_lpg_900_mj_deliv * 173,
    delta_biomass_co2e_900 = delta_biomass_mj_900_deliv * 275,
    
    delta_lpg_co2e_700 = delta_lpg_700_mj_deliv * 173,
    delta_biomass_co2e_700 = delta_biomass_mj_700_deliv * 275,
    
    delta_lpg_co2e_550 = delta_lpg_550_mj_deliv * 173,
    delta_biomass_co2e_550 = delta_biomass_mj_550_deliv * 275
  )  %>% 
  mutate(
    delta_co2e_900 = (delta_lpg_co2e_900 + delta_biomass_co2e_900)/ 1000000000000,# grams to million kg (megatonne)
    delta_co2e_700 = (delta_lpg_co2e_700 + delta_biomass_co2e_700)/ 1000000000000,
    delta_co2e_550 = (delta_lpg_co2e_550 + delta_biomass_co2e_550)/ 1000000000000
  )%>% 
  mutate(
    delta_co2e_900_usd = delta_co2e_900 * 203, 
    delta_co2e_700_usd = delta_co2e_700 * 203,
    delta_co2e_550_usd = delta_co2e_550 * 203,
    
    delta_co2e_900_usd_low = delta_co2e_900 * 100, 
    delta_co2e_700_usd_low = delta_co2e_700 * 100,
    delta_co2e_550_usd_low = delta_co2e_550 * 100,
    
    delta_co2e_900_usd_high = delta_co2e_900 * 1000, 
    delta_co2e_700_usd_high = delta_co2e_700 * 1000,
    delta_co2e_550_usd_high = delta_co2e_550 * 1000,
  ) %>% 
  mutate(
    delta_co2e_900_usd_discounted = delta_co2e_900_usd / (1 + 0.09)^(Year-2020), # 203 usd per tone in 2020, discount to 2020
    delta_co2e_700_usd_discounted = delta_co2e_700_usd / (1 + 0.09)^(Year-2020),
    delta_co2e_550_usd_discounted = delta_co2e_550_usd / (1 + 0.05)^(Year-2020),
    
    delta_co2e_900_usd_discounted_low = delta_co2e_900_usd_low / (1 + 0.09)^(Year-2020), # 203 usd per tone in 2020, discount to 2020
    delta_co2e_700_usd_discounted_low = delta_co2e_700_usd_low / (1 + 0.09)^(Year-2020),
    delta_co2e_550_usd_discounted_low = delta_co2e_550_usd_low / (1 + 0.05)^(Year-2020),
    
    delta_co2e_900_usd_discounted_high = delta_co2e_900_usd_high / (1 + 0.09)^(Year-2020), # 203 usd per tone in 2020, discount to 2020
    delta_co2e_700_usd_discounted_high = delta_co2e_700_usd_high / (1 + 0.09)^(Year-2020),
    delta_co2e_550_usd_discounted_high = delta_co2e_550_usd_high / (1 + 0.05)^(Year-2020),
  ) %>% 
  group_by(model) %>% 
  mutate(
    delta_co2e_900_usd_discounted_cumsum = cumsum(delta_co2e_900_usd_discounted),
    delta_co2e_700_usd_discounted_cumsum = cumsum(delta_co2e_700_usd_discounted),
    delta_co2e_550_usd_discounted_cumsum = cumsum(delta_co2e_550_usd_discounted),
    
    delta_co2e_900_usd_discounted_cumsum_low = cumsum(delta_co2e_900_usd_discounted_low),
    delta_co2e_700_usd_discounted_cumsum_low = cumsum(delta_co2e_700_usd_discounted_low),
    delta_co2e_550_usd_discounted_cumsum_low = cumsum(delta_co2e_550_usd_discounted_low),
    
    delta_co2e_900_usd_discounted_cumsum_high = cumsum(delta_co2e_900_usd_discounted_high),
    delta_co2e_700_usd_discounted_cumsum_high = cumsum(delta_co2e_700_usd_discounted_high),
    delta_co2e_550_usd_discounted_cumsum_high = cumsum(delta_co2e_550_usd_discounted_high),
  ) %>% 
  ungroup() %>% 
  dplyr::select(Year, 
                delta_co2e_900_usd_discounted_cumsum,
                delta_co2e_700_usd_discounted_cumsum,
                delta_co2e_550_usd_discounted_cumsum,
                
                delta_co2e_900_usd_discounted_cumsum_low,
                delta_co2e_700_usd_discounted_cumsum_low,
                delta_co2e_550_usd_discounted_cumsum_low,
                
                delta_co2e_900_usd_discounted_cumsum_high,
                delta_co2e_700_usd_discounted_cumsum_high,
                delta_co2e_550_usd_discounted_cumsum_high) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 12, 14),
    SCC = str_sub(name, nchar(name)-3,nchar(name))
  ) %>% 
  mutate(
    SCC = ifelse(SCC=="msum", "SCC\n(preferred)",
                 ifelse(SCC=="_low", "SCC\n(low)",
                        ifelse(SCC=="high", "SCC\n(high)",
                               NA))),
    price = paste0(price, " INR")
  ) %>% 
  filter(Year==2030)


india_mortality_benefits_fig <-
  ggplot() + 
  geom_errorbar(data=india_mortality_summary %>% 
                  group_by(VSL, price) %>% 
                  summarize(
                    low=quantile(value/1000000000, probs=.025),
                    high=quantile(value/1000000000, probs=.975)),
                aes(x=VSL, ymin=low, ymax=high, group=VSL), linetype="dotted", color="black", width=0.15) +
  
  geom_boxplot(data=india_mortality_summary, 
               aes(x=VSL, y=value/1000000000, group=VSL), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[1]) + 
    stat_summary(data=india_mortality_summary, 
                 aes(x=VSL, y=value/1000000000, group=VSL), fun = "mean", colour = "white", size = 2, geom = "point") +
    
    scale_x_discrete(limits=c("VSL\n(low)", "VSL\n(preferred)", "VSL\n(high)")) +
    
    ggtitle("Mortality related benefits") + 
    xlab("") + 
    ylab("USD (billions)") + 
    
    theme_bw() +
  
    facet_wrap(.~price, scales="free_y")


india_co2_benefits_fig <-
  ggplot() + 
  geom_errorbar(data=india_co2e %>% 
                  group_by(SCC, price) %>% 
                  summarize(
                    low=quantile(-value/1000, probs=.025),
                    high=quantile(-value/1000, probs=.975)),
                aes(x=SCC, ymin=low, ymax=high, group=SCC), linetype="dotted", color="black", width=0.15) +
  
  geom_boxplot(data=india_co2e, 
               aes(x=SCC, y=-value/1000, group=SCC), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[2]) + 
  stat_summary(data=india_co2e, 
               aes(x=SCC, y=-value/1000, group=SCC), fun = "mean", colour = "white", size =1, geom = "point") +
  
  scale_x_discrete(limits=c("SCC\n(low)", "SCC\n(preferred)", "SCC\n(high)")) +
  
  ggtitle("CO2e related benefits") + 
  xlab("") + 
  ylab("USD (billions)") + 
  
  theme_bw() +
  
  facet_wrap(.~price, scales="free_y")




india_monetized_sensitivity_fig <- 
  cowplot::plot_grid(
    NULL, NULL,
    india_mortality_benefits_fig,
    india_co2_benefits_fig, nrow=2, 
    rel_heights=c(0.1, 1), 
    align="hv",
    labels="Sensitivity of estimates of monetized benefits in India to higher or lower VSLs and SCCs",
    label_x = -.5
  )



# KENYA --------

kenya_mortality_summary <- 
  kenya_kg_pm_scenario_1_yr_long %>% 
  # filter(!is.na(kg_lpg_1100)) %>% 
  group_by(model) %>% 
  mutate(
    pred_1_16_deaths_difference1_VSL_discounted_cumsum = cumsum(pred_1_16_deaths_difference1_VSL_discounted),
    pred_1_8_deaths_difference1_VSL_discounted_cumsum = cumsum(pred_1_8_deaths_difference1_VSL_discounted),
    
    pred_1_16_deaths_difference1_VSL_discounted_cumsum_low = cumsum(pred_1_16_deaths_difference1_VSL_low_discounted),
    pred_1_8_deaths_difference1_VSL_discounted_cumsum_low = cumsum(pred_1_8_deaths_difference1_VSL_low_discounted),
    
    pred_1_16_deaths_difference1_VSL_discounted_cumsum_high = cumsum(pred_1_16_deaths_difference1_VSL_high_discounted),
    pred_1_8_deaths_difference1_VSL_discounted_cumsum_high = cumsum(pred_1_8_deaths_difference1_VSL_high_discounted),
  ) %>% 
  ungroup() %>% 
  filter(Year==2030) %>% 
  mutate(
    VAT0_main = pred_1_16_deaths_difference1_VSL_discounted_cumsum,
    VAT8_main = pred_1_16_deaths_difference1_VSL_discounted_cumsum - pred_1_8_deaths_difference1_VSL_discounted_cumsum,
    
    VAT0_low = pred_1_16_deaths_difference1_VSL_discounted_cumsum_low,
    VAT8_low = pred_1_16_deaths_difference1_VSL_discounted_cumsum_low - pred_1_8_deaths_difference1_VSL_discounted_cumsum_low,
    
    VAT0_high = pred_1_16_deaths_difference1_VSL_discounted_cumsum_high,
    VAT8_high = pred_1_16_deaths_difference1_VSL_discounted_cumsum_high - pred_1_8_deaths_difference1_VSL_discounted_cumsum_high,
  ) %>% 
  dplyr::select(
    Year, 
    VAT0_main, VAT8_main,
    VAT0_low, VAT8_low,
    VAT0_high, VAT8_high
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    VAT = str_sub(name, 0, 4),
    VSL = str_sub(name, 6, nchar(name))
  ) %>% 
  mutate(
    VAT = ifelse(VAT=="VAT0", "No VAT", "8% VAT"),
    VSL = ifelse(VSL=="main", "VSL\n(preferred)", 
                 ifelse(VSL=="low", "VSL\n(low)", 
                        ifelse(VSL=="high", "VSL\n(high)",NA)))
  )
  
kenya_mortality_summary

kenya_co2e <- 
  kenya_kg_pm_scenario_1_yr_long %>%  
  mutate(
    kg_lpg_consumed_observed_mj = kg_lpg_consumed_observed / 1.6 * 45,
    kg_lpg_consumed_1_16_mj = kg_lpg_consumed_1_16  / 1.6* 45,
    kg_lpg_consumed_1_8_mj = kg_lpg_consumed_1_8 / 1.6 * 45
  ) %>% 
  mutate(
    delta_lpg_mj_1_16 = kg_lpg_consumed_1_16_mj - kg_lpg_consumed_observed_mj,
    delta_lpg_mj_1_8 = kg_lpg_consumed_1_8_mj - kg_lpg_consumed_observed_mj
  ) %>% 
  mutate(
    biomass_mj_1_16 = delta_lpg_mj_1_16 * (.5/.15),
    biomass_mj_1_8 = delta_lpg_mj_1_8 * (.5/.15)
  ) %>% 
  mutate(
    delta_lpg_co2e_1_16 = delta_lpg_mj_1_16*173,
    biomass_co2e_1_16 = biomass_mj_1_16*395,
    
    delta_lpg_co2e_1_8 = delta_lpg_mj_1_8*173,
    biomass_co2e_1_8 = biomass_mj_1_8*395
  )  %>% 
  mutate(
    delta_co2e_1_16 = (delta_lpg_co2e_1_16 + biomass_co2e_1_16)/ 1000000000000,
    delta_co2e_1_8 = (delta_lpg_co2e_1_8 + biomass_co2e_1_8)/ 1000000000000
  )%>% 
  mutate(
    delta_co2e_1_16_usd = delta_co2e_1_16 * 203 * 1e6, # from Burke et al. in 2020 203 USD / ton
    delta_co2e_1_8_usd = delta_co2e_1_8 * 203 * 1e6,
    
    delta_co2e_1_16_usd_low = delta_co2e_1_16 * 100 * 1e6, 
    delta_co2e_1_8_usd_low = delta_co2e_1_8 * 100 * 1e6,
    
    delta_co2e_1_16_usd_high = delta_co2e_1_16 * 1000 * 1e6, 
    delta_co2e_1_8_usd_high = delta_co2e_1_8 * 1000 * 1e6,
  ) %>%  
  mutate(
    delta_co2e_1_16_usd_discounted = delta_co2e_1_16_usd / (1 + 0.05)^(Year-2020),
    delta_co2e_1_8_usd_discounted = delta_co2e_1_8_usd / (1 + 0.05)^(Year-2020),
    
    delta_co2e_1_16_usd_discounted_low = delta_co2e_1_16_usd_low / (1 + 0.05)^(Year-2020),
    delta_co2e_1_8_usd_discounted_low = delta_co2e_1_8_usd_low / (1 + 0.05)^(Year-2020),
    
    delta_co2e_1_16_usd_discounted_high = delta_co2e_1_16_usd_high / (1 + 0.05)^(Year-2020),
    delta_co2e_1_8_usd_discounted_high = delta_co2e_1_8_usd_high / (1 + 0.05)^(Year-2020),
  ) %>% 
  group_by(model) %>% 
  mutate(
    delta_co2e_1_16_usd_discounted_cumsum = cumsum(delta_co2e_1_16_usd_discounted),
    delta_co2e_1_8_usd_discounted_cumsum = cumsum(delta_co2e_1_8_usd_discounted),
    
    delta_co2e_1_16_usd_discounted_cumsum_low = cumsum(delta_co2e_1_16_usd_discounted_low),
    delta_co2e_1_8_usd_discounted_cumsum_low = cumsum(delta_co2e_1_8_usd_discounted_low),
    
    delta_co2e_1_16_usd_discounted_cumsum_high = cumsum(delta_co2e_1_16_usd_discounted_high),
    delta_co2e_1_8_usd_discounted_cumsum_high = cumsum(delta_co2e_1_8_usd_discounted_high),
  ) %>% 
  ungroup() %>% 
  mutate(
    VAT0_main = delta_co2e_1_16_usd_discounted_cumsum,
    VAT8_main = delta_co2e_1_16_usd_discounted_cumsum - delta_co2e_1_8_usd_discounted_cumsum,
    
    VAT0_low = delta_co2e_1_16_usd_discounted_cumsum_low,
    VAT8_low = delta_co2e_1_16_usd_discounted_cumsum_low - delta_co2e_1_8_usd_discounted_cumsum_low,
    
    VAT0_high = delta_co2e_1_16_usd_discounted_cumsum_high,
    VAT8_high = delta_co2e_1_16_usd_discounted_cumsum_high - delta_co2e_1_8_usd_discounted_cumsum_high,
  ) %>% 
  dplyr::select(Year, 
                VAT0_main,VAT8_main,
                VAT0_low,VAT8_low,
                VAT0_high,VAT8_high) %>% 
  pivot_longer(-Year) %>%
  filter(Year==2030) %>% 
  mutate(
    VAT = str_sub(name, 0, 4),
    SCC = str_sub(name, 6, nchar(name))
  ) %>% 
  mutate(
    VAT = ifelse(VAT=="VAT0", "No VAT", "8% VAT"),
    SCC = ifelse(SCC=="main", "SCC\n(preferred)", 
                 ifelse(SCC=="low", "SCC\n(low)", 
                        ifelse(SCC=="high", "SCC\n(high)",NA)))
  )

kenya_mortality_benefits_fig <-
  ggplot() + 
  geom_errorbar(data=kenya_mortality_summary %>% 
                  group_by(VSL, VAT) %>% 
                  summarize(
                    low=quantile(-value/1e9, probs=.025),
                    high=quantile(-value/1e9, probs=.975)),
                aes(x=VSL, ymin=low, ymax=high, group=VSL), linetype="dotted", color="black", width=0.15) +
  
  geom_boxplot(data=kenya_mortality_summary, 
               aes(x=VSL, y=-value/1e9, group=VSL), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[1]) + 
  stat_summary(data=kenya_mortality_summary, 
               aes(x=VSL, y=-value/1e9, group=VSL), fun = "mean", colour = "white", size = 2, geom = "point") +
  
  scale_x_discrete(limits=c("VSL\n(low)", "VSL\n(preferred)", "VSL\n(high)")) +
  
  ggtitle("Mortality related benefits") + 
  xlab("") + 
  ylab("USD (billions)") + 
  
  theme_bw() +
  
  facet_wrap(.~VAT, scales="free_y")


kenya_co2_benefits_fig <-
  ggplot() + 
  geom_errorbar(data=kenya_co2e %>% 
                  group_by(SCC, VAT) %>% 
                  summarize(
                    low=quantile(-value/1e9, probs=.025),
                    high=quantile(-value/1e9, probs=.975)),
                aes(x=SCC, ymin=low, ymax=high, group=SCC), linetype="dotted", color="black", width=0.15) +
  
  geom_boxplot(data=kenya_co2e, 
               aes(x=SCC, y=-value/1e9, group=SCC), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[2]) + 
  stat_summary(data=kenya_co2e, 
               aes(x=SCC, y=-value/1e9, group=SCC), fun = "mean", colour = "white", size = 1, geom = "point") +
  
  scale_x_discrete(limits=c("SCC\n(low)", "SCC\n(preferred)", "SCC\n(high)")) +
  
  ggtitle("CO2e related benefits") + 
  xlab("") + 
  ylab("USD (billions)") + 
  
  theme_bw() +
  
  facet_wrap(.~VAT, scales="free_y")




kenya_monetized_sensitivity_fig <- 
  cowplot::plot_grid(
    NULL, NULL,
    kenya_mortality_benefits_fig,
    kenya_co2_benefits_fig, nrow=2, 
    rel_heights=c(0.1, 1), 
    align="hv",
    labels="Sensitivity of estimates of monetized benefits in Kenya to higher or lower VSLs and SCCs",
    label_x = -.5
  )



monetized_sensitivity_combined_fig <- 
  
  plot_grid(
    ecuador_monetized_sensitivity_fig,
    india_monetized_sensitivity_fig,
    kenya_monetized_sensitivity_fig,
    nrow=3,
    align="hv"
    
  )


cowplot::ggsave2("~/Desktop/combined_inpraise_fig_sensitivity_r1.pdf",
                 monetized_sensitivity_combined_fig,
                 width = 300,
                 height = 350,
                 units = "mm",
                 dpi = 300)

# BENEFITS-COSTS RATIO FIG  --------

# Ecuador -----
ecuador_summary_costs_benefits <- 
  ecuador_summary_co2e %>% 
  dplyr::select(Year, co2e_usd_cumulative, co2e_usd_cumulative_low, co2e_usd_cumulative_high) %>% 
  left_join(ecuador_summary_costs %>% dplyr::select(Year,cost_fred_cumsum)) %>% 
  left_join(ecuador_summary_mortality %>% dplyr::select(Year, mortality_change_VSL_20yr_cumulative, mortality_change_VSL_20yr_cumulative_low, mortality_change_VSL_20yr_cumulative_high)) %>% 
  mutate(
    total_benefits = -co2e_usd_cumulative+-mortality_change_VSL_20yr_cumulative,
    total_benefits_low = -co2e_usd_cumulative_low + -mortality_change_VSL_20yr_cumulative_low,
    total_benefits_high = -co2e_usd_cumulative_high  + -mortality_change_VSL_20yr_cumulative_high
  ) %>% 
  mutate(
    benefits_vs_costs = (total_benefits*1e3)/cost_fred_cumsum,
    benefits_vs_costs_low = (total_benefits_low*1e3)/cost_fred_cumsum,
    benefits_vs_costs_high = (total_benefits_high*1e3)/cost_fred_cumsum,
  ) %>% 
  dplyr::select(Year, benefits_vs_costs, benefits_vs_costs_low, benefits_vs_costs_high) %>% 
  pivot_longer(-Year)


ecuador_benefits_vs_costs_fig <- 
  ggplot(
  
) + 
  geom_errorbar(data=ecuador_summary_costs_benefits %>% 
                  group_by(name) %>% 
                  summarize(
                    low = quantile(value, probs=.025),
                    mean = mean(value),
                    median = median(value),
                    high = quantile(value, probs=.975)
                    ), 
                aes(x=name, ymin=low, ymax=high), linetype="dotted", color="black", width=0.15) + 
  geom_boxplot(data=ecuador_summary_costs_benefits, 
               aes(x=name, y=value), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[3]) + 
  
  stat_summary(data=ecuador_summary_costs_benefits, 
               aes(x=name, y=value), fun = "mean", colour = "white", size = 2, geom = "point") +
  ylab("Ratio of benefits vs. costs") + 
  scale_x_discrete(
    limits=c("benefits_vs_costs_low", "benefits_vs_costs","benefits_vs_costs_high"),
    labels=c("Low VSL/\nLow SCC", "Preferred VSL/\nPreferred SCC","High VSL/\nHigh SCC"),
  ) + 
  xlab("") + 
  ggtitle("Ecuador") + 
  
  theme_bw()
  
# India -----


india_mortality_summary <- 
  india_kg_pm_scenario_1_yr_long %>% 
  # filter(!is.na(kg_lpg_1100)) %>% 
  group_by(model) %>% 
  mutate(
    # pred_1100_deaths_difference_VSL1 = cumsum(pred_1100_deaths_difference_VSL1),
    pred_900_deaths_difference_VSL1 = cumsum(pred_900_deaths_difference_VSL1_discounted/2.736592),
    pred_700_deaths_difference_VSL1 = cumsum(pred_700_deaths_difference_VSL1_discounted/2.736592),
    pred_550_deaths_difference_VSL1 = cumsum(pred_550_deaths_difference_VSL1_discounted/2.736592),
    
    pred_900_deaths_difference_VSL1_low = cumsum(pred_900_deaths_difference_VSL1_low_discounted/2.736592),
    pred_700_deaths_difference_VSL1_low = cumsum(pred_700_deaths_difference_VSL1_low_discounted/2.736592),
    pred_550_deaths_difference_VSL1_low = cumsum(pred_550_deaths_difference_VSL1_low_discounted/2.736592),
    
    pred_900_deaths_difference_VSL1_high = cumsum(pred_900_deaths_difference_VSL1_high_discounted/2.736592),
    pred_700_deaths_difference_VSL1_high = cumsum(pred_700_deaths_difference_VSL1_high_discounted/2.736592),
    pred_550_deaths_difference_VSL1_high = cumsum(pred_550_deaths_difference_VSL1_high_discounted/2.736592),
    
  ) %>% 
  ungroup() %>% 
  dplyr::select(Year, 
                pred_900_deaths_difference_VSL1,
                pred_700_deaths_difference_VSL1,
                pred_550_deaths_difference_VSL1,
                
                pred_900_deaths_difference_VSL1_low,
                pred_700_deaths_difference_VSL1_low,
                pred_550_deaths_difference_VSL1_low,
                
                pred_900_deaths_difference_VSL1_high,
                pred_700_deaths_difference_VSL1_high,
                pred_550_deaths_difference_VSL1_high) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 6, 8),
    VSL = str_sub(name, nchar(name)-3,nchar(name))
  ) %>% 
  mutate(
    VSL = ifelse(VSL=="VSL1", "VSL\n(preferred)",
                 ifelse(VSL=="_low", "VSL\n(low)",
                        ifelse(VSL=="high", "VSL\n(high)",
                               NA))),
    price = paste0(price, " INR")
  ) %>% 
  filter(Year==2030)


# co2e benefits -------


india_mortality_summary <- 
  india_kg_pm_scenario_1_yr_long %>% 
  # filter(!is.na(kg_lpg_1100)) %>% 
  group_by(model) %>% 
  mutate(
    # pred_1100_deaths_difference_VSL1 = cumsum(pred_1100_deaths_difference_VSL1),
    pred_900_deaths_difference_VSL1 = cumsum(pred_900_deaths_difference_VSL1_discounted/2.736592),
    pred_700_deaths_difference_VSL1 = cumsum(pred_700_deaths_difference_VSL1_discounted/2.736592),
    pred_550_deaths_difference_VSL1 = cumsum(pred_550_deaths_difference_VSL1_discounted/2.736592),
    
    pred_900_deaths_difference_VSL1_low = cumsum(pred_900_deaths_difference_VSL1_low_discounted/2.736592),
    pred_700_deaths_difference_VSL1_low = cumsum(pred_700_deaths_difference_VSL1_low_discounted/2.736592),
    pred_550_deaths_difference_VSL1_low = cumsum(pred_550_deaths_difference_VSL1_low_discounted/2.736592),
    
    pred_900_deaths_difference_VSL1_high = cumsum(pred_900_deaths_difference_VSL1_high_discounted/2.736592),
    pred_700_deaths_difference_VSL1_high = cumsum(pred_700_deaths_difference_VSL1_high_discounted/2.736592),
    pred_550_deaths_difference_VSL1_high = cumsum(pred_550_deaths_difference_VSL1_high_discounted/2.736592),
    
  ) %>% 
  ungroup() %>% 
  dplyr::select(Year, model,
                pred_900_deaths_difference_VSL1,
                pred_700_deaths_difference_VSL1,
                pred_550_deaths_difference_VSL1,
                
                pred_900_deaths_difference_VSL1_low,
                pred_700_deaths_difference_VSL1_low,
                pred_550_deaths_difference_VSL1_low,
                
                pred_900_deaths_difference_VSL1_high,
                pred_700_deaths_difference_VSL1_high,
                pred_550_deaths_difference_VSL1_high) 


india_costs_benefits <- 
  india_kg_pm_scenario_1_yr_long %>% 
  filter(!is.na(pred_1100_kg)) %>% 
  mutate(
    lpg_1100_mj = pred_1100_kg *  45, 
    lpg_900_mj = pred_900_kg * 45, 
    lpg_700_mj = pred_700_kg * 45, 
    lpg_550_mj = pred_550_kg * 45, 
  ) %>% 
  mutate(
    lpg_1100_mj_deliv = lpg_1100_mj * .5,
    lpg_900_mj_deliv = lpg_900_mj * .5,
    lpg_700_mj_deliv = lpg_700_mj * .5,
    lpg_550_mj_deliv = lpg_550_mj * .5
  ) %>% 
  mutate(
    delta_lpg_900_mj_deliv = lpg_1100_mj_deliv - lpg_900_mj_deliv,
    delta_lpg_700_mj_deliv = lpg_1100_mj_deliv - lpg_700_mj_deliv,
    delta_lpg_550_mj_deliv = lpg_1100_mj_deliv - lpg_550_mj_deliv
  ) %>% 
  mutate(
    delta_biomass_mj_900_deliv = delta_lpg_900_mj_deliv * (.5/.15),
    delta_biomass_mj_700_deliv = delta_lpg_700_mj_deliv * (.5/.15),
    delta_biomass_mj_550_deliv = delta_lpg_550_mj_deliv * (.5/.15)
  ) %>% 
  mutate(
    delta_lpg_co2e_900 = delta_lpg_900_mj_deliv * 173,
    delta_biomass_co2e_900 = delta_biomass_mj_900_deliv * 275,
    
    delta_lpg_co2e_700 = delta_lpg_700_mj_deliv * 173,
    delta_biomass_co2e_700 = delta_biomass_mj_700_deliv * 275,
    
    delta_lpg_co2e_550 = delta_lpg_550_mj_deliv * 173,
    delta_biomass_co2e_550 = delta_biomass_mj_550_deliv * 275
  )  %>% 
  mutate(
    delta_co2e_900 = (delta_lpg_co2e_900 + delta_biomass_co2e_900)/ 1000000000000,# grams to million kg (megatonne)
    delta_co2e_700 = (delta_lpg_co2e_700 + delta_biomass_co2e_700)/ 1000000000000,
    delta_co2e_550 = (delta_lpg_co2e_550 + delta_biomass_co2e_550)/ 1000000000000
  )%>% 
  mutate(
    delta_co2e_900_usd = delta_co2e_900 * 203, 
    delta_co2e_700_usd = delta_co2e_700 * 203,
    delta_co2e_550_usd = delta_co2e_550 * 203,
    
    delta_co2e_900_usd_low = delta_co2e_900 * 100, 
    delta_co2e_700_usd_low = delta_co2e_700 * 100,
    delta_co2e_550_usd_low = delta_co2e_550 * 100,
    
    delta_co2e_900_usd_high = delta_co2e_900 * 1000, 
    delta_co2e_700_usd_high = delta_co2e_700 * 1000,
    delta_co2e_550_usd_high = delta_co2e_550 * 1000,
  ) %>% 
  mutate(
    delta_co2e_900_usd_discounted = delta_co2e_900_usd / (1 + 0.09)^(Year-2020), # 203 usd per tone in 2020, discount to 2020
    delta_co2e_700_usd_discounted = delta_co2e_700_usd / (1 + 0.09)^(Year-2020),
    delta_co2e_550_usd_discounted = delta_co2e_550_usd / (1 + 0.05)^(Year-2020),
    
    delta_co2e_900_usd_discounted_low = delta_co2e_900_usd_low / (1 + 0.09)^(Year-2020), # 203 usd per tone in 2020, discount to 2020
    delta_co2e_700_usd_discounted_low = delta_co2e_700_usd_low / (1 + 0.09)^(Year-2020),
    delta_co2e_550_usd_discounted_low = delta_co2e_550_usd_low / (1 + 0.05)^(Year-2020),
    
    delta_co2e_900_usd_discounted_high = delta_co2e_900_usd_high / (1 + 0.09)^(Year-2020), # 203 usd per tone in 2020, discount to 2020
    delta_co2e_700_usd_discounted_high = delta_co2e_700_usd_high / (1 + 0.09)^(Year-2020),
    delta_co2e_550_usd_discounted_high = delta_co2e_550_usd_high / (1 + 0.05)^(Year-2020),
  ) %>% 
  mutate(
  ) %>% 
  mutate(
    diff_900 = pred_1100_kg - pred_900_kg,
    diff_700 = pred_1100_kg - pred_700_kg,
    diff_550 = pred_1100_kg - pred_550_kg
  ) %>% 
  mutate(
    subidy_900_usd_total = (.93-.76) * diff_900,
    subidy_700_usd_total = (.93-.59) * diff_900,
    subidy_550_usd_total = (.93-.47) * diff_550
  ) %>% 
  mutate(
    subidy_900_usd_total_disc = subidy_900_usd_total / (1 + 0.09)^(Year-2023),
    subidy_700_usd_total_disc = subidy_700_usd_total / (1 + 0.09)^(Year-2023),
    subidy_550_usd_total_disc = subidy_550_usd_total / (1 + 0.09)^(Year-2023)
  ) %>% 
  group_by(model) %>% 
  mutate(
    delta_co2e_900_usd_discounted_cumsum = cumsum(delta_co2e_900_usd_discounted),
    delta_co2e_700_usd_discounted_cumsum = cumsum(delta_co2e_700_usd_discounted),
    delta_co2e_550_usd_discounted_cumsum = cumsum(delta_co2e_550_usd_discounted),
    
    delta_co2e_900_usd_discounted_cumsum_low = cumsum(delta_co2e_900_usd_discounted_low),
    delta_co2e_700_usd_discounted_cumsum_low = cumsum(delta_co2e_700_usd_discounted_low),
    delta_co2e_550_usd_discounted_cumsum_low = cumsum(delta_co2e_550_usd_discounted_low),
    
    delta_co2e_900_usd_discounted_cumsum_high = cumsum(delta_co2e_900_usd_discounted_high),
    delta_co2e_700_usd_discounted_cumsum_high = cumsum(delta_co2e_700_usd_discounted_high),
    delta_co2e_550_usd_discounted_cumsum_high = cumsum(delta_co2e_550_usd_discounted_high),
  ) %>% 
  mutate(
    subidy_900_usd_total_disc_cumsum = cumsum(subidy_900_usd_total_disc),
    subidy_700_usd_total_disc_cumsum = cumsum(subidy_700_usd_total_disc),
    subidy_550_usd_total_disc_cumsum = cumsum(subidy_550_usd_total_disc)
  ) %>% 
  filter(Year==2030) %>% 
  dplyr::select(
    Year, subidy_900_usd_total_disc_cumsum, subidy_700_usd_total_disc_cumsum, subidy_550_usd_total_disc_cumsum,
    delta_co2e_900_usd_discounted, delta_co2e_700_usd_discounted, delta_co2e_550_usd_discounted,
    delta_co2e_900_usd_discounted_low, delta_co2e_700_usd_discounted_low, delta_co2e_550_usd_discounted_low,
    delta_co2e_900_usd_discounted_high, delta_co2e_700_usd_discounted_high, delta_co2e_550_usd_discounted_high,
    
  ) %>% 
  left_join(india_mortality_summary,by=c("model", "Year")) %>% 
  mutate(
    subsidy_900_ratio = (delta_co2e_900_usd_discounted + pred_900_deaths_difference_VSL1) / subidy_900_usd_total_disc_cumsum,
    subsidy_700_ratio = (delta_co2e_700_usd_discounted + pred_700_deaths_difference_VSL1) / subidy_700_usd_total_disc_cumsum,
    subsidy_550_ratio = (delta_co2e_550_usd_discounted + pred_550_deaths_difference_VSL1) / subidy_550_usd_total_disc_cumsum,
    
    subsidy_900_ratio_low = (delta_co2e_900_usd_discounted_low + pred_900_deaths_difference_VSL1_low) / subidy_900_usd_total_disc_cumsum,
    subsidy_700_ratio_low = (delta_co2e_700_usd_discounted_low + pred_700_deaths_difference_VSL1_low) / subidy_700_usd_total_disc_cumsum,
    subsidy_550_ratio_low = (delta_co2e_550_usd_discounted_low + pred_550_deaths_difference_VSL1_low) / subidy_550_usd_total_disc_cumsum,
    
    subsidy_900_ratio_high = (delta_co2e_900_usd_discounted_high + pred_900_deaths_difference_VSL1_high) / subidy_900_usd_total_disc_cumsum,
    subsidy_700_ratio_high = (delta_co2e_700_usd_discounted_high + pred_700_deaths_difference_VSL1_high) / subidy_700_usd_total_disc_cumsum,
    subsidy_550_ratio_high = (delta_co2e_550_usd_discounted_high + pred_550_deaths_difference_VSL1_high) / subidy_550_usd_total_disc_cumsum,
  ) %>% 
  ungroup() %>% 
  dplyr::select(
    Year, 
    subsidy_900_ratio, subsidy_700_ratio, subsidy_550_ratio,
    subsidy_900_ratio_low, subsidy_700_ratio_low, subsidy_550_ratio_low,
    subsidy_900_ratio_high, subsidy_700_ratio_high, subsidy_550_ratio_high
  ) %>% 
  pivot_longer(-Year) %>% 
  mutate(
    price = str_sub(name, 9, 11),
    scenario = str_sub(name, nchar(name)-4, nchar(name))
  ) %>% 
  mutate(
    price = paste0(price, " INR"),
    scenario = ifelse(scenario=="ratio", "Preferred",
                      ifelse(scenario=="o_low", "Low",
                             ifelse(scenario=="_high", "High", NA)))
  )

india_costs_benefits

india_benefits_vs_costs_fig <- 
  ggplot(
    
  ) + 
  geom_errorbar(data=india_costs_benefits %>% 
                  group_by(price, scenario) %>% 
                  summarize(
                    low = quantile(value, probs=.025),
                    mean = mean(value),
                    median = median(value),
                    high = quantile(value, probs=.975)
                  ), 
                aes(x=scenario, ymin=-low, ymax=-high, group=scenario), linetype="dotted", color="black", width=0.15) + 
  geom_boxplot(data=india_costs_benefits, 
               aes(x=scenario, y=-value, group=scenario), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[3]) + 
  
  stat_summary(data=india_costs_benefits, 
               aes(x=scenario, y=-value), fun = "mean", colour = "white", size = 2, geom = "point") +
  ylab("Ratio of benefits vs. costs") +
  scale_x_discrete(
    limits=c("Low", "Preferred","High"),
    labels=c("Low VSL/\nLow SCC", "Preferred VSL/\nPreferred SCC","High VSL/\nHigh SCC"),
  ) +
  xlab("") + 
  ggtitle("India") + 
  
  theme_bw() + 
  facet_wrap(.~price)

india_benefits_vs_costs_fig

# Kenya ------


kenya_mortality_summary <- 
  kenya_kg_pm_scenario_1_yr_long %>% 
  # filter(!is.na(kg_lpg_1100)) %>% 
  group_by(model) %>% 
  mutate(
    pred_1_16_deaths_difference1_VSL_discounted_cumsum = cumsum(pred_1_16_deaths_difference1_VSL_discounted),
    pred_1_8_deaths_difference1_VSL_discounted_cumsum = cumsum(pred_1_8_deaths_difference1_VSL_discounted),
    
    pred_1_16_deaths_difference1_VSL_discounted_cumsum_low = cumsum(pred_1_16_deaths_difference1_VSL_low_discounted),
    pred_1_8_deaths_difference1_VSL_discounted_cumsum_low = cumsum(pred_1_8_deaths_difference1_VSL_low_discounted),
    
    pred_1_16_deaths_difference1_VSL_discounted_cumsum_high = cumsum(pred_1_16_deaths_difference1_VSL_high_discounted),
    pred_1_8_deaths_difference1_VSL_discounted_cumsum_high = cumsum(pred_1_8_deaths_difference1_VSL_high_discounted),
  ) %>% 
  ungroup() %>% 
  filter(Year==2030) %>% 
  mutate(
    VAT0_main = pred_1_16_deaths_difference1_VSL_discounted_cumsum,
    VAT8_main = pred_1_16_deaths_difference1_VSL_discounted_cumsum - pred_1_8_deaths_difference1_VSL_discounted_cumsum,
    
    VAT0_low = pred_1_16_deaths_difference1_VSL_discounted_cumsum_low,
    VAT8_low = pred_1_16_deaths_difference1_VSL_discounted_cumsum_low - pred_1_8_deaths_difference1_VSL_discounted_cumsum_low,
    
    VAT0_high = pred_1_16_deaths_difference1_VSL_discounted_cumsum_high,
    VAT8_high = pred_1_16_deaths_difference1_VSL_discounted_cumsum_high - pred_1_8_deaths_difference1_VSL_discounted_cumsum_high,
  ) %>% 
  dplyr::select(
    Year, model,
    VAT0_main, VAT8_main,
    VAT0_low, VAT8_low,
    VAT0_high, VAT8_high
  ) 

kenya_mortality_summary

kenya_co2e <- 
  kenya_kg_pm_scenario_1_yr_long %>%  
  mutate(
    kg_lpg_consumed_observed_mj = kg_lpg_consumed_observed / 1.6 * 45,
    kg_lpg_consumed_1_16_mj = kg_lpg_consumed_1_16  / 1.6* 45,
    kg_lpg_consumed_1_8_mj = kg_lpg_consumed_1_8 / 1.6 * 45
  ) %>% 
  mutate(
    delta_lpg_mj_1_16 = kg_lpg_consumed_1_16_mj - kg_lpg_consumed_observed_mj,
    delta_lpg_mj_1_8 = kg_lpg_consumed_1_8_mj - kg_lpg_consumed_observed_mj
  ) %>% 
  mutate(
    biomass_mj_1_16 = delta_lpg_mj_1_16 * (.5/.15),
    biomass_mj_1_8 = delta_lpg_mj_1_8 * (.5/.15)
  ) %>% 
  mutate(
    delta_lpg_co2e_1_16 = delta_lpg_mj_1_16*173,
    biomass_co2e_1_16 = biomass_mj_1_16*395,
    
    delta_lpg_co2e_1_8 = delta_lpg_mj_1_8*173,
    biomass_co2e_1_8 = biomass_mj_1_8*395
  )  %>% 
  mutate(
    delta_co2e_1_16 = (delta_lpg_co2e_1_16 + biomass_co2e_1_16)/ 1000000000000,
    delta_co2e_1_8 = (delta_lpg_co2e_1_8 + biomass_co2e_1_8)/ 1000000000000
  )%>% 
  mutate(
    delta_co2e_1_16_usd = delta_co2e_1_16 * 203 * 1e6, # from Burke et al. in 2020 203 USD / ton
    delta_co2e_1_8_usd = delta_co2e_1_8 * 203 * 1e6,
    
    delta_co2e_1_16_usd_low = delta_co2e_1_16 * 100 * 1e6, 
    delta_co2e_1_8_usd_low = delta_co2e_1_8 * 100 * 1e6,
    
    delta_co2e_1_16_usd_high = delta_co2e_1_16 * 1000 * 1e6, 
    delta_co2e_1_8_usd_high = delta_co2e_1_8 * 1000 * 1e6,
  ) %>%  
  mutate(
    delta_co2e_1_16_usd_discounted = delta_co2e_1_16_usd / (1 + 0.05)^(Year-2020),
    delta_co2e_1_8_usd_discounted = delta_co2e_1_8_usd / (1 + 0.05)^(Year-2020),
    
    delta_co2e_1_16_usd_discounted_low = delta_co2e_1_16_usd_low / (1 + 0.05)^(Year-2020),
    delta_co2e_1_8_usd_discounted_low = delta_co2e_1_8_usd_low / (1 + 0.05)^(Year-2020),
    
    delta_co2e_1_16_usd_discounted_high = delta_co2e_1_16_usd_high / (1 + 0.05)^(Year-2020),
    delta_co2e_1_8_usd_discounted_high = delta_co2e_1_8_usd_high / (1 + 0.05)^(Year-2020),
  ) %>% 
  group_by(model) %>% 
  mutate(
    delta_co2e_1_16_usd_discounted_cumsum = cumsum(delta_co2e_1_16_usd_discounted),
    delta_co2e_1_8_usd_discounted_cumsum = cumsum(delta_co2e_1_8_usd_discounted),
    
    delta_co2e_1_16_usd_discounted_cumsum_low = cumsum(delta_co2e_1_16_usd_discounted_low),
    delta_co2e_1_8_usd_discounted_cumsum_low = cumsum(delta_co2e_1_8_usd_discounted_low),
    
    delta_co2e_1_16_usd_discounted_cumsum_high = cumsum(delta_co2e_1_16_usd_discounted_high),
    delta_co2e_1_8_usd_discounted_cumsum_high = cumsum(delta_co2e_1_8_usd_discounted_high),
  ) %>% 
  ungroup() %>% 
  mutate(
    VAT0_co2e_main = delta_co2e_1_16_usd_discounted_cumsum,
    VAT8_co2e_main = delta_co2e_1_16_usd_discounted_cumsum - delta_co2e_1_8_usd_discounted_cumsum,
    
    VAT0_co2e_low = delta_co2e_1_16_usd_discounted_cumsum_low,
    VAT8_co2e_low = delta_co2e_1_16_usd_discounted_cumsum_low - delta_co2e_1_8_usd_discounted_cumsum_low,
    
    VAT0_co2e_high = delta_co2e_1_16_usd_discounted_cumsum_high,
    VAT8_co2e_high = delta_co2e_1_16_usd_discounted_cumsum_high - delta_co2e_1_8_usd_discounted_cumsum_high,
  ) %>% 
  dplyr::select(Year, model,
                VAT0_co2e_main, VAT8_co2e_main,
                VAT0_co2e_low,VAT8_co2e_low,
                VAT0_co2e_high,VAT8_co2e_high) %>% 
  filter(Year==2030)  %>% filter(!is.na(model))%>% 
  ungroup() 

kenya_costs <- 
  kenya_kg_pm_scenario_1_yr_long %>%
  mutate(
    kg_lpg_consumed_observed =  kg_lpg_consumed_observed / 1.6,
    kg_lpg_consumed_1_16 =  kg_lpg_consumed_1_16 / 1.6,
    kg_lpg_consumed_1_8 =  kg_lpg_consumed_1_8 / 1.6
  ) %>% 
  mutate(
    diff_1_16 = kg_lpg_consumed_observed - kg_lpg_consumed_1_16,
    diff_1_8 = kg_lpg_consumed_observed - kg_lpg_consumed_1_8
  ) %>% 
  mutate(
    subidy_1_16_usd_total = (base_price*.16 -base_price) * diff_1_16,
    subidy_1_8_usd_total = (base_price*.08 -base_price) * diff_1_8
  ) %>% 
  mutate(
    subidy_1_16_usd_total_disc = subidy_1_16_usd_total / (1 + 0.09)^(Year-2023),
    subidy_1_8_usd_total_disc = subidy_1_8_usd_total / (1 + 0.09)^(Year-2023)
  ) %>% 
  group_by(model) %>% 
  mutate(
    subidy_1_16_usd_total_disc_cumsum = cumsum(subidy_1_16_usd_total_disc),
    subidy_1_8_usd_total_disc_cumsum = cumsum(subidy_1_8_usd_total_disc)
  ) %>% 
  filter(Year==2030) %>% 
  mutate(
    VAT0_cost = subidy_1_16_usd_total_disc_cumsum,
    VAT8_cost = subidy_1_16_usd_total_disc_cumsum - subidy_1_8_usd_total_disc_cumsum
  ) %>% 
  dplyr::select(
    model, Year, 
    VAT0_cost, VAT8_cost
  ) %>% 
  ungroup() %>% filter(!is.na(model))

kenya_costs_benefits <- 
  kenya_co2e %>% 
  filter(!is.na(model)) %>% 
  left_join(
      kenya_costs, by=c("model", "Year"), relationship="many-to-many"
  ) %>% 
  left_join(
    kenya_mortality_summary %>% ungroup() %>% filter(!is.na(model)), by=c("model", "Year"), relationship="many-to-many"
  ) %>% 
  mutate(
    VAT0_ratio = (VAT0_co2e_main + VAT0_main) / (VAT0_cost),
    VAT0_ratio_low = (VAT0_co2e_low + VAT0_low) / (VAT0_cost),
    VAT0_ratio_high = (VAT0_co2e_high + VAT0_high) / (VAT0_cost),
    
    VAT8_ratio = (VAT8_co2e_main + VAT8_main) / (VAT8_cost),
    VAT8_ratio_low = (VAT8_co2e_low + VAT8_low) / (VAT8_cost),
    VAT8_ratio_high = (VAT8_co2e_high + VAT8_high) / (VAT8_cost),
  ) %>% 
  dplyr::select(
    VAT0_ratio, VAT0_ratio_low, VAT0_ratio_high,
    VAT8_ratio, VAT8_ratio_low, VAT8_ratio_high
  ) %>% 
  pivot_longer(everything()) %>% 
  mutate(
    VAT = str_sub(name, 4, 4),
    scenario = str_sub(name, nchar(name)-3, nchar(name))
  ) %>% 
  mutate(
    scenario = ifelse(scenario=="atio", "Preferred",
                 ifelse(scenario=="_low", "Low",
                        ifelse(scenario=="high", "High", NA))),
    VAT = ifelse(VAT=="0", "No VAT",
                 "8% VAT")
    
  )
  

kenya_benefits_vs_costs_fig <- 
  ggplot(
    
  ) + 
  geom_errorbar(data=kenya_costs_benefits %>% 
                  group_by(VAT, scenario) %>% 
                  summarize(
                    low = quantile(value, probs=.025),
                    mean = mean(value),
                    median = median(value),
                    high = quantile(value, probs=.975)
                  ), 
                aes(x=scenario, ymin=low, ymax=high, group=scenario), linetype="dotted", color="black", width=0.15) + 
  geom_boxplot(data=kenya_costs_benefits, 
               aes(x=scenario, y=value, group=scenario), width=0.15, color="black",outlier.shape = NA, fill=met.brewer("Ingres", 3)[3]) + 
  
  stat_summary(data=kenya_costs_benefits, 
               aes(x=scenario, y=value), fun = "mean", colour = "white", size = 2, geom = "point") +
  ylab("Ratio of benefits vs. costs") +
  scale_x_discrete(
    limits=c("Low", "Preferred","High"),
    labels=c("Low VSL/\nLow SCC", "Preferred VSL/\nPreferred SCC","High VSL/\nHigh SCC"),
  ) +
  xlab("") + 
  ggtitle("Kenya") + 
  
  theme_bw() + 
  coord_cartesian(ylim=c(0, 100)) + 
  facet_wrap(.~VAT)

kenya_benefits_vs_costs_fig

benefits_costs_combined_fig <- 
  plot_grid(
    ecuador_benefits_vs_costs_fig,
    india_benefits_vs_costs_fig,
    kenya_benefits_vs_costs_fig,
    nrow=3,
    align="hv"
  )





cowplot::ggsave2("~/Desktop/combined_inpraise_fig_sensitivity_ratio.pdf",
                 benefits_costs_combined_fig,
                 width = 250,
                 height = 250,
                 units = "mm",
                 dpi = 300)
